\(~\)
\(~\)

library(Seurat)
## Attaching SeuratObject
library(grid)

so = readRDS("/mnt/market/anclab-rstudio-server/home/swil0005/projects/Wnt11/raw_data/Wnt11/wnt11-QC-clustered-sctransform.RDS")

DimPlot(so , reduction = "umap" , label = T)
grid.text("(Fig1)", 0.15, 0.85)

DimPlot(so , reduction = "umap" , group.by = "genotype")
grid.text("(Fig2)", 0.15, 0.85)

head(Idents(so) , 20)
## TGACTCCCATCATCTT AGGTTACCAGTCTGGC CTACTATCACATATGC CTCATGCTCTAACGCA 
##                0                0                0               18 
## TGAGGTTCAATAACGA GTGATGTGTAGTCCTA CACATGACATCTCAAG GTCACTCAGTTACTCG 
##               21                2                4                3 
## CACGTTCGTTCTCCCA TCATGCCCATCGCTCT CTGAGGCCAGCAGTGA TCTCTGGTCTTCCCGA 
##                1                1                1                9 
## CGACAGCGTAACCCGC CATAAGCCAGCGAACA GAGACCCAGTCGCCAC CGTAAGTGTAAGCAAT 
##                0                5                0                0 
## GAGTGTTAGCCTTCTC AAAGTCCGTAGATCCT CGAGGAACACAAGTTC GGTGTTAAGCGCACAA 
##               13                8               19                2 
## Levels: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
head(so[[]])
##                  orig.ident nCount_RNA nFeature_RNA nCount_HTO nFeature_HTO
## TGACTCCCATCATCTT wnt11mouse      22788         4784       3072            5
## AGGTTACCAGTCTGGC wnt11mouse      15595         3744       2669            5
## CTACTATCACATATGC wnt11mouse      35180         6039       5276            5
## CTCATGCTCTAACGCA wnt11mouse      20594         4508       3863            5
## TGAGGTTCAATAACGA wnt11mouse      30360         5296       3926            5
## GTGATGTGTAGTCCTA wnt11mouse      11949         3254       1532            5
##                                HTO_maxID            HTO_secondID HTO_margin
## TGACTCCCATCATCTT TS-A-53-CAGGTTGTTGTCATT TS-A-55-GTTGTGAGCACGAGA  1.7065694
## AGGTTACCAGTCTGGC TS-A-51-AACCTTTGCCACTGC TS-A-54-TATTTCCACCCGGTC  1.6277081
## CTACTATCACATATGC TS-A-54-TATTTCCACCCGGTC TS-A-52-GTCCGACTAATAGCT  2.2671312
## CTCATGCTCTAACGCA TS-A-51-AACCTTTGCCACTGC TS-A-55-GTTGTGAGCACGAGA  2.0316834
## TGAGGTTCAATAACGA TS-A-51-AACCTTTGCCACTGC TS-A-53-CAGGTTGTTGTCATT  1.8152870
## GTGATGTGTAGTCCTA TS-A-54-TATTTCCACCCGGTC TS-A-55-GTTGTGAGCACGAGA  0.8745784
##                       HTO_classification HTO_classification.global
## TGACTCCCATCATCTT TS-A-53-CAGGTTGTTGTCATT                   Singlet
## AGGTTACCAGTCTGGC TS-A-51-AACCTTTGCCACTGC                   Singlet
## CTACTATCACATATGC TS-A-54-TATTTCCACCCGGTC                   Singlet
## CTCATGCTCTAACGCA TS-A-51-AACCTTTGCCACTGC                   Singlet
## TGAGGTTCAATAACGA TS-A-51-AACCTTTGCCACTGC                   Singlet
## GTGATGTGTAGTCCTA TS-A-54-TATTTCCACCCGGTC                   Singlet
##                                  hash.ID percent.mito     S.Score  G2M.Score
## TGACTCCCATCATCTT TS-A-53-CAGGTTGTTGTCATT     1.505178 -0.04425451  0.3310395
## AGGTTACCAGTCTGGC TS-A-51-AACCTTTGCCACTGC     2.378968  0.04743516 -0.2732218
## CTACTATCACATATGC TS-A-54-TATTTCCACCCGGTC     2.151791 -0.08918952  0.3582894
## CTCATGCTCTAACGCA TS-A-51-AACCTTTGCCACTGC     3.700107 -0.37091462 -0.3498180
## TGAGGTTCAATAACGA TS-A-51-AACCTTTGCCACTGC     2.292490 -0.27135004  0.5070627
## GTGATGTGTAGTCCTA TS-A-54-TATTTCCACCCGGTC     2.426981 -0.26513688  0.3523543
##                  Phase high.mito nCount_SCT nFeature_SCT SCT_snn_res.0.2
## TGACTCCCATCATCTT   G2M     FALSE      16238         4714               2
## AGGTTACCAGTCTGGC     S     FALSE      15480         3743               0
## CTACTATCACATATGC   G2M     FALSE      16372         5006               0
## CTCATGCTCTAACGCA    G1     FALSE      16504         4448               7
## TGAGGTTCAATAACGA   G2M     FALSE      16269         4731               1
## GTGATGTGTAGTCCTA   G2M     FALSE      14626         3253               2
##                  SCT_snn_res.0.3 SCT_snn_res.0.4 SCT_snn_res.0.5
## TGACTCCCATCATCTT               1               1               1
## AGGTTACCAGTCTGGC               0               0               0
## CTACTATCACATATGC               0               0               0
## CTCATGCTCTAACGCA               8               8               8
## TGAGGTTCAATAACGA               2               2               2
## GTGATGTGTAGTCCTA               1               1               1
##                  SCT_snn_res.0.6 SCT_snn_res.0.7 SCT_snn_res.0.8
## TGACTCCCATCATCTT               0               1               0
## AGGTTACCAGTCTGGC               0               0               0
## CTACTATCACATATGC               0               0               0
## CTCATGCTCTAACGCA              12              11              18
## TGAGGTTCAATAACGA               1              19              21
## GTGATGTGTAGTCCTA               2               1               2
##                  SCT_snn_res.0.9 SCT_snn_res.1 SCT_snn_res.1.1 SCT_snn_res.1.2
## TGACTCCCATCATCTT               1             7               7               6
## AGGTTACCAGTCTGGC               0             7               7               6
## CTACTATCACATATGC               0             3               4               2
## CTCATGCTCTAACGCA              20            20              21              24
## TGAGGTTCAATAACGA              15            16              16              19
## GTGATGTGTAGTCCTA               1             0               0              10
##                  seurat_clusters  genotype
## TGACTCCCATCATCTT               6 wild-type
## AGGTTACCAGTCTGGC               6 wild-type
## CTACTATCACATATGC               2 wnt11-hom
## CTCATGCTCTAACGCA              24 wild-type
## TGAGGTTCAATAACGA              19 wild-type
## GTGATGTGTAGTCCTA              10 wnt11-hom
so@assays
## $RNA
## Assay data with 31053 features for 8065 cells
## Top 10 variable features:
##  Hbb-y, Hba-x, S100a8, Fabp4, Stfa1, Cd74, Fxyd2, S100a9, Hmox1, Ifitm1 
## 
## $HTO
## Assay data with 5 features for 8065 cells
## First 5 features:
##  TS-A-51-AACCTTTGCCACTGC, TS-A-52-GTCCGACTAATAGCT,
## TS-A-53-CAGGTTGTTGTCATT, TS-A-54-TATTTCCACCCGGTC,
## TS-A-55-GTTGTGAGCACGAGA 
## 
## $SCT
## Assay data with 18549 features for 8065 cells
## Top 10 variable features:
##  Apoe, C1qb, Hbb-bs, Hba-a1, C1qa, Fcer1g, Lyz2, Hba-a2, Pf4, Tyrobp
so@reductions
## $pca
## A dimensional reduction object with key PC_ 
##  Number of dimensions: 50 
##  Projected dimensional reduction calculated:  FALSE 
##  Jackstraw run: FALSE 
##  Computed using assay: SCT 
## 
## $umap
## A dimensional reduction object with key UMAP_ 
##  Number of dimensions: 2 
##  Projected dimensional reduction calculated:  FALSE 
##  Jackstraw run: FALSE 
##  Computed using assay: SCT
DefaultAssay(so)
## [1] "SCT"
# Number of unique groups in each meta.data column
lapply(apply(so[[]] , 2 , unique) , length)
## $orig.ident
## [1] 1
## 
## $nCount_RNA
## [1] 6897
## 
## $nFeature_RNA
## [1] 3580
## 
## $nCount_HTO
## [1] 3820
## 
## $nFeature_HTO
## [1] 1
## 
## $HTO_maxID
## [1] 5
## 
## $HTO_secondID
## [1] 5
## 
## $HTO_margin
## [1] 8060
## 
## $HTO_classification
## [1] 5
## 
## $HTO_classification.global
## [1] 1
## 
## $hash.ID
## [1] 5
## 
## $percent.mito
## [1] 8039
## 
## $S.Score
## [1] 8065
## 
## $G2M.Score
## [1] 8060
## 
## $Phase
## [1] 3
## 
## $high.mito
## [1] 1
## 
## $nCount_SCT
## [1] 2142
## 
## $nFeature_SCT
## [1] 2478
## 
## $SCT_snn_res.0.2
## [1] 12
## 
## $SCT_snn_res.0.3
## [1] 15
## 
## $SCT_snn_res.0.4
## [1] 15
## 
## $SCT_snn_res.0.5
## [1] 16
## 
## $SCT_snn_res.0.6
## [1] 19
## 
## $SCT_snn_res.0.7
## [1] 21
## 
## $SCT_snn_res.0.8
## [1] 22
## 
## $SCT_snn_res.0.9
## [1] 23
## 
## $SCT_snn_res.1
## [1] 25
## 
## $SCT_snn_res.1.1
## [1] 25
## 
## $SCT_snn_res.1.2
## [1] 28
## 
## $seurat_clusters
## [1] 28
## 
## $genotype
## [1] 2
which(unlist(lapply(apply(so[[]] , 2 , unique) , length)) == 22)
## SCT_snn_res.0.8 
##              25
unique(so[[]][ , 25])
##  [1] 0  18 21 2  4  3  1  9  5  13 8  19 7  12 11 6  20 10 16 14 17 15
## Levels: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
so.sub = subset(so , idents = c("1" , "5" , "13" , "15" , "18" , "21"))
so.sub
## An object of class Seurat 
## 49607 features across 1895 samples within 3 assays 
## Active assay: SCT (18549 features, 3000 variable features)
##  2 other assays present: RNA, HTO
##  2 dimensional reductions calculated: pca, umap
DimPlot(so.sub , reduction = "umap" , label = T)
grid.text("(Fig3)", 0.15, 0.85)

DimPlot(so.sub , reduction = "umap" , group.by = "genotype")
grid.text("(Fig4)", 0.15, 0.85)

DimPlot(so.sub , reduction = "umap" , group.by = "Phase")
grid.text("(Fig5)", 0.15, 0.85)

\(~\)
\(~\)

plot = DimPlot(so.sub , reduction = "umap" , label = T)
selected = CellSelector(plot = plot)
selected

so.small = CellSelector(plot = plot , object = so.sub , ident = 'SelectedCells')


head(colnames(so.small[[]]))

DimPlot(so.small)

unique(so.small[[]][,25])

lapply(lapply(apply(so.small[[]] , 2 , unique) , function(x) x %in% "SelectedCells") , which)

Idents(so.sub)

Idents(so.small)

levels(so.sub)

levels(droplevels(Idents(so.sub)))

levels(so.small)

levels(droplevels(Idents(so.small)))




so.small = subset(so.small , idents = "SelectedCells")

DimPlot(so.small)

Idents(so.small) = "SCT_snn_res.0.8"

DimPlot(so.small , label = T , label.size = 6)

\(~\)
\(~\)

library(Seurat)
library(SingleCellExperiment)
library(slingshot)
library(RColorBrewer)
library(bioc2020trajectories)
## Warning: replacing previous import 'HDF5Array::path' by 'igraph::path' when
## loading 'bioc2020trajectories'
## Warning: replacing previous import 'igraph::compose' by 'purrr::compose' when
## loading 'bioc2020trajectories'
## Warning: replacing previous import 'igraph::simplify' by 'purrr::simplify' when
## loading 'bioc2020trajectories'
library(ggbeeswarm)
library(ggthemes)
library(BiocParallel)
library(tradeSeq)
library(future)
library(tidyverse)


so.small = readRDS("/mnt/market/anclab-rstudio-server/home/mpir0002/Slingshot/so.small.rds")

\(~\)
\(~\)

Trajectories on Clusters 1, 5, 13, 15, 18, and 21

DimPlot(so.small , reduction = "umap" , label = T , label.size = 6)
grid.text("(Fig6)", 0.15, 0.85)

sce <- as.SingleCellExperiment(so.small)
sce
## class: SingleCellExperiment 
## dim: 18549 1847 
## metadata(0):
## assays(2): counts logcounts
## rownames(18549): Xkr4 Sox17 ... CAAA01118383.1 CAAA01147332.1
## rowData names(6): sct.detection_rate sct.gmean ...
##   sct.residual_variance sct.variable
## colnames(1847): CTCATGCTCTAACGCA TGAGGTTCAATAACGA ... ATTCACTGTCCAGAAG
##   AGACAAAAGTACTGTC
## colData names(32): orig.ident nCount_RNA ... genotype ident
## reducedDimNames(2): PCA UMAP
## altExpNames(0):
head(Idents(so.small))
## CTCATGCTCTAACGCA TGAGGTTCAATAACGA CACGTTCGTTCTCCCA TCATGCCCATCGCTCT 
##               18               21                1                1 
## CTGAGGCCAGCAGTGA CATAAGCCAGCGAACA 
##                1                5 
## Levels: 1 5 13 15 18 21
slingshot.1to18 <- slingshot(sce, reducedDim = 'UMAP', clusterLabels = sce$SCT_snn_res.0.8 , start.clus = "1" , end.clus = "18")
## Using full covariance matrix
slingshot.1to18
## class: SingleCellExperiment 
## dim: 18549 1847 
## metadata(0):
## assays(2): counts logcounts
## rownames(18549): Xkr4 Sox17 ... CAAA01118383.1 CAAA01147332.1
## rowData names(6): sct.detection_rate sct.gmean ...
##   sct.residual_variance sct.variable
## colnames(1847): CTCATGCTCTAACGCA TGAGGTTCAATAACGA ... ATTCACTGTCCAGAAG
##   AGACAAAAGTACTGTC
## colData names(35): orig.ident nCount_RNA ... slingPseudotime_1
##   slingPseudotime_2
## reducedDimNames(2): PCA UMAP
## altExpNames(0):
sds.1to18 <- SlingshotDataSet(slingshot.1to18)
sds.1to18
## class: SlingshotDataSet 
## 
##  Samples Dimensions
##     1847          2
## 
## lineages: 2 
## Lineage1: 1  21  5  13  18  
## Lineage2: 1  21  5  13  15  
## 
## curves: 2 
## Curve1: Length: 10.129   Samples: 1672.77
## Curve2: Length: 10.861   Samples: 1710.37
slingshot.1to15 <- slingshot(sce, reducedDim = 'UMAP', clusterLabels = sce$SCT_snn_res.0.8 , start.clus = "1" , end.clus = "15")
## Using full covariance matrix
slingshot.1to15
## class: SingleCellExperiment 
## dim: 18549 1847 
## metadata(0):
## assays(2): counts logcounts
## rownames(18549): Xkr4 Sox17 ... CAAA01118383.1 CAAA01147332.1
## rowData names(6): sct.detection_rate sct.gmean ...
##   sct.residual_variance sct.variable
## colnames(1847): CTCATGCTCTAACGCA TGAGGTTCAATAACGA ... ATTCACTGTCCAGAAG
##   AGACAAAAGTACTGTC
## colData names(35): orig.ident nCount_RNA ... slingPseudotime_1
##   slingPseudotime_2
## reducedDimNames(2): PCA UMAP
## altExpNames(0):
sds.1to15 <- SlingshotDataSet(slingshot.1to15)
sds.1to15
## class: SlingshotDataSet 
## 
##  Samples Dimensions
##     1847          2
## 
## lineages: 2 
## Lineage1: 1  21  5  13  18  
## Lineage2: 1  21  5  13  15  
## 
## curves: 2 
## Curve1: Length: 10.129   Samples: 1672.77
## Curve2: Length: 10.861   Samples: 1710.37

\(~\)
\(~\)

Plotting the trajectories

colors <- rainbow(100, alpha = 1)

plot(reducedDims(slingshot.1to18)$UMAP, 
     col = colors[cut(slingshot.1to18$slingPseudotime_1,breaks=100)],
     pch=16, asp = 1, main = "colored_basedon_Pseudotime1_cluster1to18")
lines(SlingshotDataSet(slingshot.1to18), lwd=2)
grid.text("(Fig7)", 0.15, 0.85)

plot(reducedDims(slingshot.1to18)$UMAP, 
     col = colors[cut(slingshot.1to18$slingPseudotime_2,breaks=100)],
     pch=16, asp = 1, main = "colored_basedon_Pseudotime2_cluster1to18")
lines(SlingshotDataSet(slingshot.1to18), lwd=2)
grid.text("(Fig8)", 0.15, 0.85)

# Show the cells coloured by clusters
plot(reducedDims(slingshot.1to18)$UMAP, col = brewer.pal(9,'Set1')[droplevels(slingshot.1to18$SCT_snn_res.0.8)], pch=16, 
     main = "colored_basedon_Clusters_cluster1to18")
lines(SlingshotDataSet(slingshot.1to18), lwd=2)
grid.text("(Fig9)", 0.15, 0.85)

\(~\)
\(~\)

colors <- rainbow(100, alpha = 1)

plot(reducedDims(slingshot.1to15)$UMAP, 
     col = colors[cut(slingshot.1to15$slingPseudotime_1,breaks=100)],
     pch=16, asp = 1, main = "colored_basedon_Pseudotime1_cluster1to15")
lines(SlingshotDataSet(slingshot.1to15), lwd=2)
grid.text("(Fig10)", 0.15, 0.85)

plot(reducedDims(slingshot.1to15)$UMAP, 
     col = colors[cut(slingshot.1to15$slingPseudotime_2,breaks=100)],
     pch=16, asp = 1, main = "colored_basedon_Pseudotime2_cluster1to15")
lines(SlingshotDataSet(slingshot.1to15), lwd=2)
grid.text("(Fig11)", 0.15, 0.85)

# Show the cells coloured by clusters
plot(reducedDims(slingshot.1to15)$UMAP, col = brewer.pal(9,'Set1')[droplevels(slingshot.1to15$SCT_snn_res.0.8)], pch=16, 
     main = "colored_basedon_Clusters_cluster1to15")
lines(SlingshotDataSet(slingshot.1to15), lwd=2)
grid.text("(Fig12)", 0.15, 0.85)

\(~\)
\(~\)
\(~\)
\(~\)

Trajectories on Clusters 1, 5, 13, and 21

so.small2 = subset(so.small , idents = c("1" , "5" , "13" , "21"))
DimPlot(so.small2 , reduction = "umap" , label = T , label.size = 6) 
grid.text("(Fig13)", 0.15, 0.85)

sce2 = as.SingleCellExperiment(so.small2)

slingshot.1to13 = slingshot(sce2 , reducedDim = 'UMAP', clusterLabels = sce2$SCT_snn_res.0.8 , start.clus = "1" , end.clus = "13")
## Using full covariance matrix
slingshot.1to13
## class: SingleCellExperiment 
## dim: 18549 1650 
## metadata(0):
## assays(2): counts logcounts
## rownames(18549): Xkr4 Sox17 ... CAAA01118383.1 CAAA01147332.1
## rowData names(6): sct.detection_rate sct.gmean ...
##   sct.residual_variance sct.variable
## colnames(1650): TGAGGTTCAATAACGA CACGTTCGTTCTCCCA ... ATTCACTGTCCAGAAG
##   AGACAAAAGTACTGTC
## colData names(34): orig.ident nCount_RNA ... slingClusters
##   slingPseudotime_1
## reducedDimNames(2): PCA UMAP
## altExpNames(0):
sds.1to13 = SlingshotDataSet(slingshot.1to13)
sds.1to13
## class: SlingshotDataSet 
## 
##  Samples Dimensions
##     1650          2
## 
## lineages: 1 
## Lineage1: 1  21  5  13  
## 
## curves: 1 
## Curve1: Length: 9.1458   Samples: 1650

\(~\)
\(~\)

Plotting the trajectories

colors <- rainbow(100, alpha = 1)

plot(reducedDims(slingshot.1to13)$UMAP, 
     col = colors[cut(slingshot.1to13$slingPseudotime_1,breaks=100)],
     pch=16, asp = 1, main = "colored_basedon_Pseudotime1_cluster1to13")
lines(SlingshotDataSet(slingshot.1to13), lwd=2)
grid.text("(Fig14)", 0.15, 0.75)

# Show the cells coloured by clusters
plot(reducedDims(slingshot.1to13)$UMAP, col = brewer.pal(9,'Set1')[droplevels(slingshot.1to13$SCT_snn_res.0.8)], pch=16, 
     main = "colored_basedon_Clusters_cluster1to13")
lines(SlingshotDataSet(slingshot.1to13), lwd=2)
grid.text("(Fig15)", 0.15, 0.75)

\(~\)
\(~\)

Compute an imbalance score to show whether nearby cells have the same condition of not

scores = imbalance_score(rd = reducedDims(sce)$UMAP, cl = colData(sce)$genotype,
        k = 20, smooth = 40)

imbalance_data = merge(as.data.frame(reducedDims(sce)$UMAP), as.data.frame(scores), by = 0, all = TRUE)
ggplot() + geom_point(aes(x = UMAP_1,y = UMAP_2,colour = scaled_scores),data=imbalance_data, size = 3) + scale_colour_viridis_c(option = "plasma") + theme_linedraw() + ggtitle("Imbalance of Wild-type and Wnt11KO cells across UMAP space")
grid.text("(Fig16)", 0.15, 0.85)

\(~\)
\(~\)

suppressPackageStartupMessages({
        library(scales)
        library(viridis); library(UpSetR)
        library(pheatmap); library(msigdbr)
        library(fgsea); library(knitr)
        library(ggplot2); library(gridExtra)
        library(tradeSeq)
})

\(~\)
\(~\)

Distribution of cells based on cell-cyle and genotype

shuffle <- sample(ncol(sce))
layout(matrix(1:2, nrow = 1))
par(mar = c(4.5,4,1,1))

plot(reducedDims(sce)$UMAP[shuffle, ],
     asp = 1, pch = 16, xlab = "UMAP-1", ylab = "UMAP-2",
     col = alpha(c(1:2)[factor(colData(sce)$genotype)][shuffle], alpha = .5))
legend("topright", pch = 16, col = 1:2, bty = "n", 
       legend = levels(factor(colData(sce)$genotype)))
grid.text("(Fig17)", 0.15, 0.75)

plot(reducedDims(sce)$UMAP[shuffle, ], asp = 1, pch = 16, xlab = "UMAP-1", ylab = "UMAP-2", 
     col = alpha(c(3,4,5)[factor(colData(sce)$Phase)][shuffle], alpha = .5))
legend("topright", pch = 16, col = 3:5, bty = "n", legend = levels(factor(colData(sce)$Phase)))
grid.text("(Fig18)", 0.25, 0.75)

\(~\)
\(~\)

Plot the distribution of cells per cluster across pseudotimes.

slingshot.1to18$SCT_snn_res.0.8 = droplevels(slingshot.1to18$SCT_snn_res.0.8)
slingshot.1to15$SCT_snn_res.0.8 = droplevels(slingshot.1to15$SCT_snn_res.0.8)
slingshot.1to13$SCT_snn_res.0.8 = droplevels(slingshot.1to13$SCT_snn_res.0.8)


ggplot(as.data.frame(colData(slingshot.1to18)),
       aes(x = slingshot.1to18$slingPseudotime_1,
           y = forcats::fct_relevel(slingshot.1to18$SCT_snn_res.0.8,"1", "5", "13", "15", "18" , "21"), colour = SCT_snn_res.0.8)) +
        geom_quasirandom(groupOnX = F) +
        scale_color_tableau() + theme_classic() +
        xlab("Slingshot pseudotime") + ylab("Clusters") +
        ggtitle("Aggregated") + NoLegend()
## Warning: Removed 111 rows containing missing values (position_quasirandom).
grid.text("(Fig19)", 0.15, 0.85)

summary(slingshot.1to18$slingPseudotime_1)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   0.000   1.210   3.002   3.540   5.427  10.129     111
so.wild = subset(so.small , subset = genotype == "wild-type")
so.wild
## An object of class Seurat 
## 49607 features across 827 samples within 3 assays 
## Active assay: SCT (18549 features, 3000 variable features)
##  2 other assays present: RNA, HTO
##  2 dimensional reductions calculated: pca, umap
so.wnt11 = subset(so.small , subset = genotype == "wnt11-hom")
so.wnt11
## An object of class Seurat 
## 49607 features across 1020 samples within 3 assays 
## Active assay: SCT (18549 features, 3000 variable features)
##  2 other assays present: RNA, HTO
##  2 dimensional reductions calculated: pca, umap
sce.wild = as.SingleCellExperiment(so.wild)
sce.wnt11 = as.SingleCellExperiment(so.wnt11)

slingshot.wild.1to18 <- slingshot(sce.wild, reducedDim = 'UMAP', clusterLabels = sce.wild$SCT_snn_res.0.8 , start.clus = "1" , end.clus = "18")
## Using full covariance matrix
slingshot.wild.1to18
## class: SingleCellExperiment 
## dim: 18549 827 
## metadata(0):
## assays(2): counts logcounts
## rownames(18549): Xkr4 Sox17 ... CAAA01118383.1 CAAA01147332.1
## rowData names(6): sct.detection_rate sct.gmean ...
##   sct.residual_variance sct.variable
## colnames(827): CTCATGCTCTAACGCA TGAGGTTCAATAACGA ... GGGCTCAAGCATCTTG
##   CGAATTGGTGGATCGA
## colData names(35): orig.ident nCount_RNA ... slingPseudotime_1
##   slingPseudotime_2
## reducedDimNames(2): PCA UMAP
## altExpNames(0):
slingshot.wnt11.1to18 <- slingshot(sce.wnt11, reducedDim = 'UMAP', clusterLabels = sce.wnt11$SCT_snn_res.0.8 , start.clus = "1" , end.clus = "18")
## Using full covariance matrix
slingshot.wnt11.1to18
## class: SingleCellExperiment 
## dim: 18549 1020 
## metadata(0):
## assays(2): counts logcounts
## rownames(18549): Xkr4 Sox17 ... CAAA01118383.1 CAAA01147332.1
## rowData names(6): sct.detection_rate sct.gmean ...
##   sct.residual_variance sct.variable
## colnames(1020): CATAAGCCAGCGAACA GAGTGTTAGCCTTCTC ... ATTCACTGTCCAGAAG
##   AGACAAAAGTACTGTC
## colData names(35): orig.ident nCount_RNA ... slingPseudotime_1
##   slingPseudotime_2
## reducedDimNames(2): PCA UMAP
## altExpNames(0):
ggplot(as.data.frame(colData(slingshot.wild.1to18)),
       aes(x = slingshot.wild.1to18$slingPseudotime_1,
           y = forcats::fct_relevel(slingshot.wild.1to18$SCT_snn_res.0.8,"1", "5", "13", "15", "18" , "21"), colour = SCT_snn_res.0.8)) +
        geom_quasirandom(groupOnX = F) +
        scale_color_tableau() + theme_classic() +
        xlab("Pseudotime1") + ylab("Clusters") +
        ggtitle("Wild Type") + NoLegend()
## Warning: Removed 35 rows containing missing values (position_quasirandom).
grid.text("(Fig20)", 0.15, 0.85)

ggplot(as.data.frame(colData(slingshot.wnt11.1to18)),
       aes(x = slingshot.wnt11.1to18$slingPseudotime_1,
           y = forcats::fct_relevel(slingshot.wnt11.1to18$SCT_snn_res.0.8,"1", "5", "13", "15", "18" , "21"), colour = SCT_snn_res.0.8)) +
        geom_quasirandom(groupOnX = F) +
        scale_color_tableau() + theme_classic() +
        xlab("Pseudotime1") + ylab("Clusters") +
        ggtitle("Wnt11KO") + NoLegend()
## Warning: Removed 80 rows containing missing values (position_quasirandom).
grid.text("(Fig21)", 0.15, 0.85)

ggplot(as.data.frame(colData(slingshot.1to18)),
       aes(x = slingshot.1to18$slingPseudotime_1,
           y = forcats::fct_relevel(slingshot.1to18$SCT_snn_res.0.8,
            "1", "5", "13", "15", "18" , "21"), 
           colour = SCT_snn_res.0.8)) +
        geom_quasirandom(groupOnX = F) +
        facet_grid(rows = vars(genotype)) +
        scale_color_tableau() + theme_classic() +
        xlab("Slingshot pseudotime") + ylab("Clusters") +
        ggtitle("") + NoLegend()
## Warning: Removed 30 rows containing missing values (position_quasirandom).
## Warning: Removed 81 rows containing missing values (position_quasirandom).
grid.text("(Fig22)", 0.15, 0.85)

\(~\)
\(~\)

Statistical test between distribution of pseudotime1 and 2 in wild-type and wnt11

# On pseudotime one
ks.test(slingPseudotime(slingshot.wild.1to18)[, 1], slingPseudotime(slingshot.wnt11.1to18)[, 1])
## Warning in ks.test(slingPseudotime(slingshot.wild.1to18)[, 1],
## slingPseudotime(slingshot.wnt11.1to18)[, : p-value will be approximate in the
## presence of ties
## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  slingPseudotime(slingshot.wild.1to18)[, 1] and slingPseudotime(slingshot.wnt11.1to18)[, 1]
## D = 0.10589, p-value = 0.0001303
## alternative hypothesis: two-sided
# On pseudotime two
ks.test(slingPseudotime(slingshot.wild.1to18)[, 2], slingPseudotime(slingshot.wnt11.1to18)[, 2])
## Warning in ks.test(slingPseudotime(slingshot.wild.1to18)[, 2],
## slingPseudotime(slingshot.wnt11.1to18)[, : p-value will be approximate in the
## presence of ties
## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  slingPseudotime(slingshot.wild.1to18)[, 2] and slingPseudotime(slingshot.wnt11.1to18)[, 2]
## D = 0.099284, p-value = 0.000381
## alternative hypothesis: two-sided

\(~\)
\(~\)

Boxplots

summary(slingPseudotime(slingshot.wild.1to18))[,1]
##                                                                             
## "Min.   : 0.000  " "1st Qu.: 1.019  " "Median : 2.841  " "Mean   : 3.275  " 
##                                                          
## "3rd Qu.: 5.228  " "Max.   :10.365  "     "NA's   :35  "
summary(slingPseudotime(slingshot.wnt11.1to18))[,1]
##                                                                             
## "Min.   : 0.000  " "1st Qu.: 1.525  " "Median : 3.376  " "Mean   : 3.949  " 
##                                                          
## "3rd Qu.: 6.012  " "Max.   :10.199  "     "NA's   :80  "
boxplot(list(wild = slingPseudotime(slingshot.wild.1to18)[,1] , wnt11 = slingPseudotime(slingshot.wnt11.1to18)[,1]) , 
        horizontal = T,  
        notch  = T,   # Confidence interval for median
        main   = "Pseudotime One",
        sub    = "",
        xlab   = "",
        col    = "gray" )
grid.text("(Fig23)", 0.15, 0.85)

boxplot(list(wild = slingPseudotime(slingshot.wild.1to18)[,2] , wnt11 = slingPseudotime(slingshot.wnt11.1to18)[,2]) , 
        horizontal = T,  
        notch  = T,   # Confidence interval for median
        main   = "Pseudotime two",
        sub    = "",
        xlab   = "",
        col    = "gray50" )
grid.text("(Fig24)", 0.15, 0.85)

\(~\)
\(~\)

QQplots

QQplot shows that there is a inharmony between the wild-type and wnt11 cells in the later pseudotimes.

par(mfrow = c(1,2))

qqplot(slingPseudotime(slingshot.wild.1to18)[,1] , slingPseudotime(slingshot.wnt11.1to18)[,1] , 
       main = "Quantile-Quantile Plot Psudotime One",
       xlab = "wild-type Percentiles",
       ylab = "Wnt11 Percentiles",
       col.axis = 'blue', col.lab = 'red', cex.axis = 1, cex.lab = 1 , cex.main = 1)
grid.text("(Fig25)", 0.15, 0.85)

qqplot(slingPseudotime(slingshot.wild.1to18)[,2] , slingPseudotime(slingshot.wnt11.1to18)[,2] , 
       main = "Quantile-Quantile Plot Psudotime two",
       xlab = "wild-type Percentiles",
       ylab = "Wnt11 Percentiles",
       col.axis = 'blue', col.lab = 'red', cex.axis = 1, cex.lab = 1 , cex.main = 1)
grid.text("(Fig26)", 0.25, 0.85)

dev.off()
## null device 
##           1

\(~\)
\(~\)

Density Plots

dens = density(as.data.frame(na.omit(slingPseudotime(slingshot.wild.1to18)))[,1])
plot(dens, frame =  F , col = "steelblue", main = "Wild-type Psudotime One") 
polygon(dens, col = "steelblue")
grid.text("(Fig27)", 0.15, 0.75)

dens = density(as.data.frame(na.omit(slingPseudotime(slingshot.wnt11.1to18)))[,1])
plot(dens, frame = FALSE, col = "indianred1", main = "Wnt11 Psudotime One") 
polygon(dens, col = "indianred1")
grid.text("(Fig28)", 0.15, 0.75)

dens = density(as.data.frame(na.omit(slingPseudotime(slingshot.wild.1to18)))[,1])
plot(dens, frame =  F , col = "steelblue", main = "Psudotime One" , xlab = "") 
grid.text("(Fig29)", 0.15, 0.75)

dens = density(as.data.frame(na.omit(slingPseudotime(slingshot.wnt11.1to18)))[,1])
lines(dens, col = "indianred1", main = "")

\(~\)
\(~\)

dens = density(as.data.frame(na.omit(slingPseudotime(slingshot.wild.1to18)))[,2])
plot(dens, frame =  F , col = "steelblue", main = "Wild-type Psudotime Two") 
polygon(dens, col = "steelblue")
grid.text("(Fig30)", 0.15, 0.75)

dens = density(as.data.frame(na.omit(slingPseudotime(slingshot.wnt11.1to18)))[,2])
plot(dens, frame = FALSE, col = "indianred1", main = "Wnt11 Psudotime Two") 
polygon(dens, col = "indianred1")
grid.text("(Fig31)", 0.15, 0.75)

dens = density(as.data.frame(na.omit(slingPseudotime(slingshot.wild.1to18)))[,2])
plot(dens, frame =  F , col = "steelblue", main = "Psudotime Two" , xlab = "") 
grid.text("(Fig32)", 0.15, 0.75)

dens = density(as.data.frame(na.omit(slingPseudotime(slingshot.wnt11.1to18)))[,2])
lines(dens, col = "indianred1", main = "")

\(~\)
\(~\)

Differential Expression Analysis

For each condition (wild and wnt11), a smooth average expression profile along pseudotime will be estimated for each gene, using a negative binomial generalized additive model (NB-GAM). But first, we need to Evaluate the optimal number of knots required for fitGAM.

plan("multiprocess", workers = 18)
## Warning in supportsMulticoreAndRStudio(...): [ONE-TIME WARNING] Forked
## processing ('multicore') is not supported when running R from RStudio
## because it is considered unstable. For more details, how to control forked
## processing or not, and how to silence this warning in future R sessions, see ?
## parallelly::supportsMulticore
BPPARAMS <- BiocParallel::bpparam()
BPPARAMS$workers <- 12
icMat = evaluateK(slingshot.1to18, conditions = factor(colData(slingshot.1to18)$genotype, 
   ordered = FALSE), nGenes = 500, k = 3:10, parallel=TRUE, BPPARAM = BPPARAMS)


Figure33

\(~\)

In the left panel, boxplots, the skewness in the distribution of gene-level AIC values changes as the number of knots increases from eight knots. In addition, the Average AIC levels off as it reaches the eighth knot.

\(~\)
\(~\)

GAM.1to18 <- fitGAM(counts = slingshot.1to18@assays@data$counts %>% as.matrix(), sds = sds.1to18, conditions = factor(colData(slingshot.1to18)$genotype, ordered = FALSE), nknots = 8)

\(~\)
\(~\)

Perform statistical test to check whether gene expression is constant across pseudotime within a lineage

To assess significant changes in gene expression as a function of pseudotime within each lineage, we use the “associationTest” (Wald test), which tests the null hypothesis that gene expression is not a function of pseudotime, i.e., whether the estimated smoothers are significantly varying as a function of pseudotime within each lineage. The lineages=TRUE argument specifies that we would like the results for each lineage separately, asides from the default global test, which tests for significant associations across all lineages in the trajectory simultaneously. Further, we specify a log2 fold change cut-off to test against using the l2fc argument.

rowData(GAM.1to18)$assocRes <- associationTest(GAM.1to18, lineages = TRUE, l2fc = log2(2))
assocRes <- rowData(GAM.1to18)$assocRes

all(assocRes$meanLogFC > 0)
## [1] TRUE

\(~\)
\(~\)

Selecting genes with fdr-adjusted p-value less than five percent on lineage1 for both conditions.

wildGenes <-  rownames(assocRes)[
        which(p.adjust(assocRes$`pvalue_lineage1_conditionwild-type`, "fdr") <= 0.05)
]
wnt11Genes <-  rownames(assocRes)[
        which(p.adjust(assocRes$`pvalue_lineage1_conditionwnt11-hom`, "fdr") <= 0.05)
]

length(wildGenes) 
## [1] 391
length(wnt11Genes) 
## [1] 358

\(~\)
\(~\)

Genes that are common between wild and wnt11

intersect(wildGenes , wnt11Genes)
##   [1] "Gsta3"         "Col3a1"        "Col5a2"        "Fn1"          
##   [5] "Ugt1a9"        "Arl4c"         "Tnfrsf11a"     "Tmem37"       
##   [9] "Plekha6"       "Cfh"           "4930500M09Rik" "Cenpf"        
##  [13] "Irf6"          "Vim"           "Pax8"          "Gsn"          
##  [17] "Zeb2"          "Dpp4"          "Ttn"           "Serping1"     
##  [21] "Rapsn"         "Gm10804"       "Rcn1"          "Prnp"         
##  [25] "Myl9"          "Mafb"          "Ube2c"         "Bgn"          
##  [29] "Cited1"        "Maged2"        "Tmsb4x"        "Ccna2"        
##  [33] "Gucy1a1"       "Fga"           "S100a6"        "S100a10"      
##  [37] "Ecm1"          "Hist2h2ac"     "Pdzk1"         "Olfml3"       
##  [41] "Arhgap29"      "Abca4"         "Npnt"          "Penk"         
##  [45] "Tpm2"          "Gm12678"       "Slc5a9"        "Cdc20"        
##  [49] "Cdca8"         "Eva1b"         "Gm12999"       "Sfn"          
##  [53] "Aunip"         "Mfap2"         "Chd5"          "Mmp23"        
##  [57] "Ttll10"        "Cenpa"         "Emilin1"       "Ldb2"         
##  [61] "Cckar"         "Pdgfra"        "Igfbp7"        "Gm33050"      
##  [65] "Bmp3"          "Plac8"         "Arhgap24"      "Pcolce"       
##  [69] "Rasl11a"       "Col1a2"        "Cald1"         "Tmem140"      
##  [73] "Rarres2"       "Npy"           "Gprin3"        "Ndnf"         
##  [77] "Asprv1"        "Bcam"          "Cblc"          "Prodh2"       
##  [81] "Fxyd7"         "Klk10"         "Emp3"          "Fah"          
##  [85] "AI314278"      "Prss23"        "Stard10"       "Rassf10"      
##  [89] "Acsm1"         "Crym"          "Tgfb1i1"       "H19"          
##  [93] "Cdkn1c"        "Ccnd1"         "Akap12"        "Gm47889"      
##  [97] "Upb1"          "Cnn2"          "Timp3"         "Srgap1"       
## [101] "Eef1akmt3"     "Hmgb2"         "BC030500"      "Sh2d4a"       
## [105] "Pde4c"         "Elmo3"         "Clec18a"       "Snai3"        
## [109] "D830030K20Rik" "Bmp4"          "Clu"           "Pdlim2"       
## [113] "Izumo1r"       "Ntm"           "Tagln"         "Hspb2"        
## [117] "Cryab"         "Slc35f2"       "Islr2"         "Pclaf"        
## [121] "Tpm1"          "Anxa2"         "Dag1"          "Selenom"      
## [125] "4921536K21Rik" "Emid1"         "Col23a1"       "Pmp22"        
## [129] "Tmem102"       "Pimreg"        "Hic1"          "Lhx1"         
## [133] "Lhx1os"        "Car4"          "Col1a1"        "Top2a"        
## [137] "Kcnj16"        "Birc5"         "Hist1h1b"      "Hist1h2ap"    
## [141] "Hist1h2ae"     "Hist1h4d"      "Hist1h1e"      "Hist1h2ab"    
## [145] "Hist1h1a"      "Gmnn"          "Gm17750"       "Bhmt2"        
## [149] "Ccnb1"         "Rrm2"          "Scin"          "Mis18bp1"     
## [153] "Meg3"          "Agxt2"         "Npr3"          "Osr2"         
## [157] "Cthrc1"        "Lgals1"        "Igfbp6"        "Cldn5"        
## [161] "Stfa3"         "Fstl1"         "Pcp4"          "Ripk4"        
## [165] "Clic5"         "Uhrf1"         "Lbh"           "Cdc42ep3"     
## [169] "Six2"          "Dsg2"          "Mep1b"         "Ccdc178"      
## [173] "Cd248"         "Slc22a19"      "Cpn1"          "Emx2"

\(~\)
\(~\)

Visualization of number of DE genes that present in either both or one cell type.

UpSetR::upset(fromList(list(wild = wildGenes, wnt11 = wnt11Genes)))
grid.text("(Fig34)", 0.15, 0.85)

\(~\)
\(~\)

DEGs that present only in wild-type

wildGenes[!(wildGenes %in% wnt11Genes)]
##   [1] "Gm16894"       "8430432A02Rik" "Plcd4"         "Inpp5d"       
##   [5] "Gm26720"       "Cdh20"         "Rassf5"        "Gm29629"      
##   [9] "Fmod"          "Mybph"         "Myog"          "Ppfia4"       
##  [13] "Rgs16"         "Pbx1"          "Slamf1"        "Tlr5"         
##  [17] "Prox1os"       "Cd46"          "Acbd7"         "Cacna1b"      
##  [21] "Tor4a"         "Tubb4b"        "Gm996"         "Lcn2"         
##  [25] "Ifih1"         "Spc25"         "Rapgef4"       "Fibin"        
##  [29] "Rhov"          "Nusap1"        "Gm14133"       "Tpx2"         
##  [33] "Hnf4aos"       "Aurka"         "Gpc3"          "Xlr4b"        
##  [37] "Trpc5"         "Lhfpl1"        "Frmpd4"        "Fabp4"        
##  [41] "Gm5103"        "Mme"           "Golim4"        "Gm15998"      
##  [45] "Hfe2"          "5330417C22Rik" "Fam166b"       "Stra6l"       
##  [49] "Smc2"          "Adamtsl1"      "Kif2c"         "Ccdc24"       
##  [53] "Zc3h12a"       "Matn1"         "Ncmap"         "Sh2d5"        
##  [57] "2700016F22Rik" "Slc25a34"      "Gpr157"        "Gm26840"      
##  [61] "Tmem88b"       "Lrrd1"         "Gm8773"        "Sema3a"       
##  [65] "Fgfbp1"        "Nipal1"        "Ptpn13"        "Abcg3"        
##  [69] "Idua"          "Plcxd1"        "Selplg"        "Mlxipl"       
##  [73] "Fzd9"          "Gm15498"       "Gper1"         "Mmd2"         
##  [77] "2900089D17Rik" "Irf5"          "Rncr4"         "Trh"          
##  [81] "Frmd4b"        "Syn2"          "Cxcl12"        "D330020A13Rik"
##  [85] "Gm38910"       "Apold1"        "Bhlhe41"       "Pnmal1"       
##  [89] "Cd177"         "Rinl"          "Gm26810"       "Nphs1"        
##  [93] "Cd33"          "Gm45669"       "Kcnj14"        "1700011D18Rik"
##  [97] "Prc1"          "Cemip"         "Me3"           "Hbb-bt"       
## [101] "Hbb-bs"        "Gm26690"       "D7Ertd443e"    "Gm44732"      
## [105] "Rgs17"         "Adgrg6"        "Cenpw"         "Mfsd4b3"      
## [109] "Tspan15"       "Cdk1"          "Plk5"          "Gna15"        
## [113] "Gm32687"       "Eid3"          "Ntn4"          "Myl6"         
## [117] "Efnb2"         "Proscos"       "Mfap3l"        "Gm15991"      
## [121] "Lzts1"         "Iqcm"          "Gm10645"       "Gm45767"      
## [125] "Gm3636"        "Gm3739"        "Duxbl3"        "Stab1"        
## [129] "Eddm3b"        "Rnase6"        "Trac"          "Cma1"         
## [133] "Loxl2"         "Itm2b"         "Igsf9b"        "Gm15520"      
## [137] "Adamts15"      "Slc37a2"       "Gm10658"       "9230112J17Rik"
## [141] "Unc13c"        "Gcm1"          "Fam46a"        "Fbxw15"       
## [145] "Dlec1"         "Ccdc13"        "Cdcp1"         "Plek"         
## [149] "Hba-a1"        "Hba-a2"        "Hmmr"          "Gabrb2"       
## [153] "Sox30"         "Flt4"          "Hs3st3b1"      "Aurkb"        
## [157] "Tm4sf5"        "Rilp"          "Ccl12"         "Wfdc17"       
## [161] "Sept4"         "Dgkeos"        "Sp6"           "Csf3"         
## [165] "Gast"          "Icam2"         "Pecam1"        "Abca6"        
## [169] "Fn3krp"        "Akr1c21"       "Hist1h3h"      "Hist1h2bc"    
## [173] "Agtr1a"        "Irf4"          "Gm26514"       "A330048O09Rik"
## [177] "Pfn3"          "Gm48488"       "Ctsl"          "Gm47123"      
## [181] "BC048507"      "Dhfr"          "Tmem171"       "Gm10734"      
## [185] "Fam228a"       "Srp54c"        "Tdrd9"         "A530016L24Rik"
## [189] "Rgs22"         "Pick1.1"       "Gm16059"       "1810041L15Rik"
## [193] "Gm26798"       "Shank3"        "Pde1b"         "Ciita"        
## [197] "Aifm3"         "Rtl10"         "Il1rap"        "Epha3"        
## [201] "Cldn6"         "Prss27"        "Slc9a3r2"      "Itpr3"        
## [205] "Scube3"        "Platr17"       "Trem2"         "Apobec2"      
## [209] "Arhgap28"      "Mkx"           "Rnf125"        "Cd14"         
## [213] "Gnal"          "Rin1"          "Trpm3"

\(~\)

DEGs that present only in wnt11

wnt11Genes[!(wnt11Genes %in% wildGenes)]
##   [1] "Eya1"          "D630023F18Rik" "Pid1"          "Col6a3"       
##   [5] "Cdk18"         "Tmem81"        "Angptl1"       "Atp1b1"       
##   [9] "Gm37422"       "Rd3"           "Nrarp"         "Phyhd1"       
##  [13] "Pip5kl1"       "Lzts3"         "Plcb1"         "Jag1"         
##  [17] "Hnf4a"         "5031425F14Rik" "Cybb"          "Rtl8a"        
##  [21] "Rtl8b"         "Gm14827"       "Gpm6b"         "Cldn11"       
##  [25] "Anxa5"         "Tm4sf1"        "Enpep"         "Lef1"         
##  [29] "Arid3c"        "Col15a1"       "Acnat1"        "Aldob"        
##  [33] "Kif12"         "Nfib"          "Ccdc17"        "Kcnq4"        
##  [37] "Tmem54"        "Gm12976"       "Fam229a"       "Wnt4"         
##  [41] "Hes5"          "Tmem52"        "Gm26894"       "Gm43660"      
##  [45] "Tacc3"         "Gm15477"       "Gabra2"        "Afm"          
##  [49] "Cdkl2"         "Tmem150c"      "Gm26711"       "Oasl2"        
##  [53] "Wnt16"         "Lrrc4"         "Podxl"         "Svopl"        
##  [57] "Epha1"         "Aoc1"          "Tacstd2"       "Eva1a"        
##  [61] "Mitf"          "Pparg"         "Clec4a2"       "Slco1a6"      
##  [65] "Cdc42ep5"      "Phldb3"        "Aspdh"         "Slc17a7"      
##  [69] "Chrna7"        "Agbl1"         "Wdr93"         "Gm44949"      
##  [73] "Tm6sf1"        "Slco2b1"       "Hbb-y"         "Spon1"        
##  [77] "Atp2a1"        "Aldoa.1"       "Gm44759"       "Ifitm3"       
##  [81] "Gm16982"       "Cd81"          "Lama2"         "Lama4"        
##  [85] "Lilrb4a"       "Srgn"          "Col6a2"        "Dcn"          
##  [89] "Atp2b1"        "Alx1"          "Lrriq1"        "Nav3"         
##  [93] "Ptprr"         "Gm36041"       "Mcemp1"        "Tgfbr3l"      
##  [97] "Col4a1"        "Proz"          "Gm32050"       "Enpp6"        
## [101] "Lrrc25"        "Il15"          "Nfix"          "Cdh16"        
## [105] "Gm33023"       "Slc35f3"       "Cfap70"        "Gm30108"      
## [109] "Wnt5a"         "Vstm4"         "Arhgap22"      "Nrg3"         
## [113] "Fermt2"        "Gm34934"       "1700011H14Rik" "Cebpe"        
## [117] "Myh7"          "Dct"           "Gria4"         "Nxpe2"        
## [121] "Arhgap20os"    "Crabp1"        "Islr"          "Trpc1"        
## [125] "Als2cl"        "Hba-x"         "Ebf1"          "Acsl6"        
## [129] "Anxa6"         "Slc5a10"       "Tnk1"          "Slc6a4"       
## [133] "Hnf1b"         "Hoxb1"         "Nags"          "Cygb"         
## [137] "Gm11728"       "CT010465.1"    "Nid1"          "Hist1h1d"     
## [141] "Hist1h3e"      "Hist1h3c"      "Dcdc2a"        "Dsp"          
## [145] "Rnf144b"       "Id4"           "Cks2"          "Cdhr2"        
## [149] "Gas1"          "Zfp455"        "Foxd1"         "Gpx8"         
## [153] "Arl4a"         "Rhoj"          "Gm26571"       "Fbln5"        
## [157] "AC163040.1"    "Dlk1"          "Amn"           "Rai14"        
## [161] "Pdzd2"         "Snhg18"        "Ncald"         "Galr3"        
## [165] "B130046B21Rik" "Smagp"         "Cdc45"         "Nrros"        
## [169] "Hcls1"         "Sytl3"         "Dll1"          "Tmem8"        
## [173] "Grm4"          "Gm16279"       "Ttbk1"         "1700061G19Rik"
## [177] "Pdgfrb"        "Ablim3"        "Myo5b"         "Acta2"        
## [181] "Ifit3"         "Scd3"

\(~\)
\(~\)

Heatmap of wild DE genes based on mean smoother. It Gets smoothers estimated by tradeSeq along a grid. This function does not return fitted values but rather the predicted mean smoother, for a user-defined grid of points. If tidy = FALSE and the trajectory consists of 2 lineages and nPoints=100, then the returned matrix will have 2*100 columns, where the first 100 correspond to the first lineage and columns 101-200 to the second lineage.

yhatSmooth <- predictSmooth(GAM.1to18, gene = wildGenes, nPoints = 100, tidy = FALSE)
heatSmooth <- pheatmap(t(scale(t(yhatSmooth[, 1:100]))),
                       cluster_cols = F, show_rownames = T, 
                       show_colnames = F, fontsize_row = 2 , main = "Wild-type DEGs on lineage1")
grid.text("(Fig35)", 0.15, 0.85)

\(~\)
\(~\)

Dendrogram of predicted smooth average expression

c1 <- sort(cutree(heatSmooth$tree_row, k = 6))
c1
##         Gsta3        Inpp5d         Arl4c     Tnfrsf11a        Tmem37 
##             1             1             1             1             1 
##        Rassf5       Plekha6          Fmod         Mybph          Myog 
##             1             1             1             1             1 
##         Rgs16        Slamf1          Tlr5          Irf6         Acbd7 
##             1             1             1             1             1 
##           Gsn          Dpp4         Ifih1           Ttn         Rapsn 
##             1             1             1             1             1 
##          Rcn1          Rhov          Prnp       Gm14133         Xlr4b 
##             1             1             1             1             1 
##         Fabp4           Mme 5330417C22Rik         Abca4          Npnt 
##             1             1             1             1             1 
##       Fam166b        Stra6l        Slc5a9        Ccdc24       Gm12999 
##             1             1             1             1             1 
##           Sfn         Ncmap 2700016F22Rik        Ttll10        Gm8773 
##             1             1             1             1             1 
##        Fgfbp1         Cckar        Nipal1        Igfbp7      Arhgap24 
##             1             1             1             1             1 
##        Ptpn13         Abcg3          Idua        Plcxd1        Mlxipl 
##             1             1             1             1             1 
##       Gm15498         Gper1          Irf5         Rncr4       Tmem140 
##             1             1             1             1             1 
##        Frmd4b          Syn2 D330020A13Rik        Apold1          Cblc 
##             1             1             1             1             1 
##       Gm26810         Klk10           Fah        Prss23           Me3 
##             1             1             1             1             1 
##       Stard10        Hbb-bt        Hbb-bs       Gm26690    D7Ertd443e 
##             1             1             1             1             1 
##       Gm44732         Rgs17        Adgrg6       Gm47889          Upb1 
##             1             1             1             1             1 
##          Plk5     Eef1akmt3       Proscos      BC030500        Mfap3l 
##             1             1             1             1             1 
##       Gm15991         Pde4c          Iqcm         Elmo3       Clec18a 
##             1             1             1             1             1 
##         Snai3 D830030K20Rik        Gm3636        Gm3739         Stab1 
##             1             1             1             1             1 
##          Bmp4        Eddm3b        Rnase6          Trac          Cma1 
##             1             1             1             1             1 
##           Clu         Itm2b       Izumo1r      Adamts15       Slc37a2 
##             1             1             1             1             1 
##         Hspb2       Slc35f2       Gm10658        Unc13c          Dag1 
##             1             1             1             1             1 
##        Fbxw15         Dlec1         Cdcp1         Emid1        Hba-a1 
##             1             1             1             1             1 
##        Hba-a2        Gabrb2      Hs3st3b1        Tm4sf5          Car4 
##             1             1             1             1             1 
##         Sept4           Sp6         Abca6        Kcnj16        Fn3krp 
##             1             1             1             1             1 
##      Hist1h3h        Agtr1a A330048O09Rik          Ctsl      BC048507 
##             1             1             1             1             1 
##       Gm17750       Tmem171       Gm10734          Scin          Osr2 
##             1             1             1             1             1 
##         Rgs22        Cthrc1 1810041L15Rik         Pde1b         Aifm3 
##             1             1             1             1             1 
##        Il1rap         Ripk4         Cldn6        Prss27         Itpr3 
##             1             1             1             1             1 
##       Platr17         Trem2       Apobec2           Mkx          Dsg2 
##             1             1             1             1             1 
##       Ccdc178          Cd14         Trpm3          Cpn1          Emx2 
##             1             1             1             1             1 
##       Gm16894 8430432A02Rik         Plcd4       Gm26720         Cdh20 
##             2             2             2             2             2 
##        Ppfia4 4930500M09Rik       Prox1os          Cd46       Cacna1b 
##             2             2             2             2             2 
##         Tor4a         Gm996          Lcn2       Rapgef4        Cited1 
##             2             2             2             2             2 
##         Trpc5        Lhfpl1        Frmpd4        Gm5103       Gm15998 
##             2             2             2             2             2 
##          Hfe2          Tpm2       Zc3h12a         Matn1         Sh2d5 
##             2             2             2             2             2 
##         Mfap2      Slc25a34        Gpr157          Chd5       Gm26840 
##             2             2             2             2             2 
##       Tmem88b         Lrrd1        Selplg          Fzd9          Mmd2 
##             2             2             2             2             2 
## 2900089D17Rik        Gprin3        Asprv1       Bhlhe41        Pnmal1 
##             2             2             2             2             2 
##         Cd177          Rinl       Gm45669        Kcnj14          Emp3 
##             2             2             2             2             2 
## 1700011D18Rik         Cemip          Crym         Gna15       Gm32687 
##             2             2             2             2             2 
##          Eid3         Lzts1       Gm10645        Duxbl3        Igsf9b 
##             2             2             2             2             2 
##         Islr2 9230112J17Rik        Fam46a        Ccdc13 4921536K21Rik 
##             2             2             2             2             2 
##          Plek         Sox30          Flt4          Rilp         Ccl12 
##             2             2             2             2             2 
##        Wfdc17        Dgkeos          Csf3          Gast         Icam2 
##             2             2             2             2             2 
##          Irf4       Gm26514          Pfn3       Gm48488       Gm47123 
##             2             2             2             2             2 
##       Fam228a        Srp54c         Tdrd9 A530016L24Rik       Gm16059 
##             2             2             2             2             2 
##       Gm26798        Shank3         Ciita         Stfa3          Six2 
##             2             2             2             2             2 
##        Rnf125          Gnal          Rin1        Col3a1        Col5a2 
##             2             2             2             3             3 
##           Fn1           Cfh          Pbx1          Zeb2      Serping1 
##             3             3             3             3             3 
##         Fibin          Myl9          Gpc3           Bgn        Maged2 
##             3             3             3             3             3 
##        Tmsb4x       Gucy1a1        S100a6       S100a10          Ecm1 
##             3             3             3             3             3 
##        Olfml3          Penk      Adamtsl1         Eva1b        Sema3a 
##             3             3             3             3             3 
##       Emilin1          Ldb2        Pdgfra         Plac8        Pcolce 
##             3             3             3             3             3 
##        Col1a2         Cald1       Rarres2         Fxyd7       Tgfb1i1 
##             3             3             3             3             3 
##        Cdkn1c          Cnn2        Pdlim2           Ntm         Anxa2 
##             3             3             3             3             3 
##       Selenom       Col23a1          Hic1        Col1a1          Meg3 
##             3             3             3             3             3 
##        Lgals1         Epha3           Lbh      Cdc42ep3         Cd248 
##             3             3             3             3             3 
##        Ugt1a9       Gm29629       Gm10804       Hnf4aos           Fga 
##             4             4             4             4             4 
##         Pdzk1       Gm12678       Gm33050           Trh       Gm38910 
##             4             4             4             4             4 
##        Prodh2          Cd33      AI314278         Acsm1       Mfsd4b3 
##             4             4             4             4             4 
##       Gm45767       Gm15520          Gcm1       Tmem102        Pecam1 
##             4             4             4             4             4 
##       Akr1c21     Hist1h2bc         Bhmt2         Agxt2       Pick1.1 
##             4             4             4             4             4 
##         Rtl10         Mep1b      Slc22a19         Cenpf        Tubb4b 
##             4             4             4             5             5 
##         Spc25        Nusap1          Tpx2         Ube2c         Aurka 
##             5             5             5             5             5 
##         Ccna2     Hist2h2ac          Smc2         Kif2c         Cdc20 
##             5             5             5             5             5 
##         Cdca8         Aunip         Cenpa          Prc1         Cenpw 
##             5             5             5             5             5 
##          Cdk1         Hmgb2         Pclaf          Hmmr         Aurkb 
##             5             5             5             5             5 
##        Pimreg         Top2a         Birc5      Hist1h1b     Hist1h2ap 
##             5             5             5             5             5 
##     Hist1h2ae      Hist1h4d      Hist1h1e     Hist1h2ab      Hist1h1a 
##             5             5             5             5             5 
##          Gmnn          Dhfr         Ccnb1          Rrm2      Mis18bp1 
##             5             5             5             5             5 
##         Uhrf1           Vim          Pax8          Mafb        Golim4 
##             5             6             6             6             6 
##      Arhgap29         Mmp23          Bmp3       Rasl11a           Npy 
##             6             6             6             6             6 
##          Ndnf        Cxcl12          Bcam         Nphs1       Rassf10 
##             6             6             6             6             6 
##           H19         Ccnd1        Akap12       Tspan15         Timp3 
##             6             6             6             6             6 
##          Ntn4        Srgap1          Myl6         Efnb2        Sh2d4a 
##             6             6             6             6             6 
##         Loxl2         Tagln         Cryab          Tpm1         Pmp22 
##             6             6             6             6             6 
##          Lhx1        Lhx1os          Npr3        Igfbp6         Cldn5 
##             6             6             6             6             6 
##         Fstl1          Pcp4      Slc9a3r2        Scube3         Clic5 
##             6             6             6             6             6 
##      Arhgap28 
##             6
table(c1)
## c1
##   1   2   3   4   5   6 
## 150  88  47  28  38  40
plot(heatSmooth$tree_row , sub = "" , main = "" , xlab = "" , cex = 0.3)
grid.text("(Fig36)", 0.15, 0.85)

\(~\)
\(~\)

Visualizing the heatmap for the actual fitted values in sorted cells based on pseudotimes. In other words, cells with larger pseudotime are toward the right-side of the heatmap. In addition, Clustering were applied only on rows (genes).

conditions <- GAM.1to18$tradeSeq$conditions

pt1 <- slingshot.1to18$slingPseudotime_1

## Based on fitted values

yhatCell <- predictCells(GAM.1to18, gene=wildGenes)

dim(yhatCell)
## [1]  391 1847
yhatCellwild <- yhatCell[ , conditions == "wild-type"]

yhatCellwild = yhatCellwild[ , !is.na(pt1[conditions == "wild-type"])]

dim(yhatCellwild)
## [1] 391 797
pt1.wild = pt1[conditions == "wild-type"][!is.na(pt1[conditions == "wild-type"])]


# Order according to pseudotime
# Remove NA pseudotimes
# Separate cells that are in pseudotime zero from the other cells. Then, combines them together to have zero-pseudotime cells at the begining of the matrix, then sorted-non-zero-psuedotime cells. 

sum(pt1.wild == 0)
## [1] 72
yhatCellwildzero = yhatCellwild[ , pt1.wild == 0]

yhatCellwild.nonzero = yhatCellwild[ , pt1.wild != 0]

dim(yhatCellwild.nonzero)
## [1] 391 725
dim(yhatCellwildzero)
## [1] 391  72
dim(yhatCellwild)
## [1] 391 797
pt1.wild.zero = pt1.wild[pt1.wild == 0]

pt1.wild.nonzero = pt1.wild[pt1.wild != 0]

pt1.wild.nonzero.sorted = sort(pt1.wild.nonzero, decreasing=FALSE)

yhatCellwild.nonzero.sorted = yhatCellwild.nonzero[ , match(pt1.wild.nonzero.sorted, pt1.wild.nonzero)]

mat = cbind(yhatCellwildzero , yhatCellwild.nonzero.sorted)

dim(mat)
## [1] 391 797
pheatmap(t(scale(t(mat))), cluster_cols = FALSE,
         show_rownames = T, show_colnames=FALSE , fontsize_row = 2 , main = "Wild-type DEGs across lineage1, wild-type sorted cells based on lineage 1")
grid.text("(Fig37)", 0.15, 0.85)

\(~\)
\(~\)

GSEA based on wild-type pseudotimes on lineage1

## C5 category is according to gene ontology grouping: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4707969/pdf/nihms-743907.pdf
geneSets <- msigdbr(species = "Mus musculus", category = "C5", subcategory = "BP")
### filter background to only include genes that we assessed.

geneSets <- geneSets[geneSets$gene_symbol %in% names(GAM.1to18),]
m_list <- geneSets %>% split(x = .$gene_symbol, f = .$gs_name)

stats <- assocRes$`waldStat_lineage1_conditionwild-type`
names(stats) <- rownames(assocRes)
eaRes <- fgsea(pathways = m_list, stats = stats, nperm = 5e4, minSize = 10)
## Warning in fgsea(pathways = m_list, stats = stats, nperm = 50000, minSize = 10):
## You are trying to run fgseaSimple. It is recommended to use fgseaMultilevel. To
## run fgseaMultilevel, you need to remove the nperm argument in the fgsea function
## call.
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (8.43% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
## Warning in fgseaSimple(...): There were 26 pathways for which P-values were not
## calculated properly due to unbalanced gene-level statistic values
eaRes = eaRes %>% arrange(pval)

head(eaRes[,c(1:3,8)] , 40)
##                                                                                                  pathway
##  1: GOBP_NEGATIVE_REGULATION_OF_VIRAL_INDUCED_CYTOPLASMIC_PATTERN_RECOGNITION_RECEPTOR_SIGNALING_PATHWAY
##  2:                                                                         GOBP_RIG_I_SIGNALING_PATHWAY
##  3:                                                           GOBP_REGULATION_OF_RIG_I_SIGNALING_PATHWAY
##  4:          GOBP_REGULATION_OF_VIRAL_INDUCED_CYTOPLASMIC_PATTERN_RECOGNITION_RECEPTOR_SIGNALING_PATHWAY
##  5:                                                GOBP_NEGATIVE_REGULATION_OF_DEFENSE_RESPONSE_TO_VIRUS
##  6:                            GOBP_NEGATIVE_REGULATION_OF_TRANSCRIPTION_BY_COMPETITIVE_PROMOTER_BINDING
##  7:                                                                             GOBP_HISTONE_METHYLATION
##  8:                                          GOBP_MYD88_INDEPENDENT_TOLL_LIKE_RECEPTOR_SIGNALING_PATHWAY
##  9:                                                          GOBP_RESPONSE_TO_LEUKEMIA_INHIBITORY_FACTOR
## 10:                                                          GOBP_POSITIVE_REGULATION_OF_RRNA_PROCESSING
## 11:                                                                      GOBP_NON_MOTILE_CILIUM_ASSEMBLY
## 12:                                                  GOBP_NEGATIVE_REGULATION_OF_MYOTUBE_DIFFERENTIATION
## 13:                 GOBP_POSITIVE_REGULATION_OF_TRANSCRIPTION_ELONGATION_FROM_RNA_POLYMERASE_II_PROMOTER
## 14:                                                          GOBP_TOLL_LIKE_RECEPTOR_4_SIGNALING_PATHWAY
## 15:                                                         GOBP_NEGATIVE_REGULATION_OF_DEFENSE_RESPONSE
## 16:                                                            GOBP_CYTOPLASMIC_MICROTUBULE_ORGANIZATION
## 17:                                                          GOBP_CEREBRAL_CORTEX_NEURON_DIFFERENTIATION
## 18:                                                                      GOBP_HISTONE_H3_K27_METHYLATION
## 19:                                              GOBP_NEGATIVE_REGULATION_OF_RESPONSE_TO_BIOTIC_STIMULUS
## 20:                                                        GOBP_INOSITOL_TRISPHOSPHATE_METABOLIC_PROCESS
## 21:                                                   GOBP_LIPOPOLYSACCHARIDE_MEDIATED_SIGNALING_PATHWAY
## 22:                                                                   GOBP_REGULATION_OF_RRNA_PROCESSING
## 23:                                       GOBP_POSITIVE_REGULATION_OF_SMOOTH_MUSCLE_CELL_DIFFERENTIATION
## 24:                                               GOBP_POSITIVE_REGULATION_OF_NEUROTRANSMITTER_SECRETION
## 25:                                              GOBP_POSITIVE_REGULATION_OF_INTERFERON_GAMMA_PRODUCTION
## 26:                                                                   GOBP_NUCLEAR_MEMBRANE_ORGANIZATION
## 27:                                                         GOBP_CIRCADIAN_REGULATION_OF_GENE_EXPRESSION
## 28:                                                                       GOBP_RESPONSE_TO_MAGNESIUM_ION
## 29:          GOBP_REGULATION_OF_NUCLEAR_TRANSCRIBED_MRNA_CATABOLIC_PROCESS_DEADENYLATION_DEPENDENT_DECAY
## 30:                            GOBP_REGULATION_OF_VASCULAR_ASSOCIATED_SMOOTH_MUSCLE_CELL_DIFFERENTIATION
## 31:                             GOBP_POSITIVE_REGULATION_OF_G_PROTEIN_COUPLED_RECEPTOR_SIGNALING_PATHWAY
## 32:                                                GOBP_POSITIVE_REGULATION_OF_SMOOTH_MUSCLE_CONTRACTION
## 33:                                                     GOBP_REGULATION_OF_LONG_TERM_SYNAPTIC_DEPRESSION
## 34:                                            GOBP_MYD88_DEPENDENT_TOLL_LIKE_RECEPTOR_SIGNALING_PATHWAY
## 35:                                GOBP_POSITIVE_REGULATION_OF_SIGNAL_TRANSDUCTION_BY_P53_CLASS_MEDIATOR
## 36:                                               GOBP_POSITIVE_REGULATION_OF_NEUROTRANSMITTER_TRANSPORT
## 37:                                GOBP_POSITIVE_REGULATION_OF_RELEASE_OF_CYTOCHROME_C_FROM_MITOCHONDRIA
## 38:                                       GOBP_POSITIVE_REGULATION_OF_ENDOTHELIAL_CELL_APOPTOTIC_PROCESS
## 39:                                             GOBP_POSITIVE_REGULATION_OF_ALCOHOL_BIOSYNTHETIC_PROCESS
## 40:                                  GOBP_POSITIVE_REGULATION_OF_CELLULAR_CARBOHYDRATE_METABOLIC_PROCESS
##                                                                                                  pathway
##            pval      padj                                 leadingEdge
##  1: 0.001138057 0.8789942                                 Rnf125,Tkfc
##  2: 0.001641642 0.8789942                                Rnf125,Ddx60
##  3: 0.001749201 0.8789942                                Rnf125,Ddx60
##  4: 0.001804258 0.8789942                                Rnf125,Ddx60
##  5: 0.002103703 0.8789942                                Rnf125,Mill2
##  6: 0.003080147 0.8789942                                Bhlhe41,Muc1
##  7: 0.003379932 0.8789942                                  Chd5,Smyd1
##  8: 0.003560071 0.8789942                                   Cd14,Cd40
##  9: 0.003639927 0.8789942                                Rnf125,Tex14
## 10: 0.003952569 0.8789942   Trmt112,Wdr43,Riok2,Dimt1,Utp15,Riok1,...
## 11: 0.004399912 0.8789942                               Ccdc13,Cep350
## 12: 0.004756417 0.8789942                              Bhlhe41,Notch1
## 13: 0.005376344 0.8789942     Tcerg1,Gtf2f1,Brd4,Eapp,Cdk13,Cdk12,...
## 14: 0.005500110 0.8789942                                    Cd14,Lbp
## 15: 0.005719886 0.8789942                       Rnf125,Gper1,Serping1
## 16: 0.006439871 0.8789942                                 Ccdc13,Mapt
## 17: 0.006504629 0.8789942                                   Chd5,Emx1
## 18: 0.006674917 0.8789942                                  Chd5,Phf19
## 19: 0.006779864 0.8789942                             Rnf125,Serping1
## 20: 0.008088004 0.8789942                                 Gper1,P2ry6
## 21: 0.008119838 0.8789942                                 Cd14,Sigirr
## 22: 0.008415147 0.8789942   Trmt112,Wdr43,Kat2b,Riok2,Dimt1,Utp15,...
## 23: 0.008423258 0.8789942                                   Gper1,Eng
## 24: 0.010326848 0.8789942                                 Gper1,Nlgn1
## 25: 0.010379792 0.8789942                                  Cd14,Wnt5a
## 26: 0.010495205 0.8789942                                 Gper1,Reep4
## 27: 0.010919782 0.8789942                                 Bhlhe41,Id4
## 28: 0.011137163 0.8789942                                  Cd14,Ccnd1
## 29: 0.011220196 0.8789942 Nanos1,Cpeb3,Zfp36l2,Tnrc6b,Ago2,Pabpc1,...
## 30: 0.011320984 0.8789942                                   Gper1,Eng
## 31: 0.011399964 0.8789942                                 Gper1,Tmod2
## 32: 0.011525642 0.8789942                                  Gper1,Oxtr
## 33: 0.012215820 0.8789942                                 Shank3,Ager
## 34: 0.012280491 0.8789942                              Cd14,Tlr5,Cd36
## 35: 0.012481968 0.8789942                                   Chd5,Hic1
## 36: 0.013531174 0.8789942                                 Gper1,Itgb1
## 37: 0.013602448 0.8789942                                   Gper1,Hrk
## 38: 0.014905591 0.8789942                                 Gper1,Cd248
## 39: 0.015216866 0.8789942                                  Gper1,Wnt4
## 40: 0.015219696 0.8789942                                 Gper1,Actn3
##            pval      padj                                 leadingEdge

\(~\)
\(~\)

GSEA based on wnt11 pseudotimes on lineage1

stats <- assocRes$`waldStat_lineage1_conditionwnt11-hom`
names(stats) <- rownames(assocRes)
eaRes <- fgsea(pathways = m_list, stats = stats, nperm = 5e4, minSize = 10)
## Warning in fgsea(pathways = m_list, stats = stats, nperm = 50000, minSize = 10):
## You are trying to run fgseaSimple. It is recommended to use fgseaMultilevel. To
## run fgseaMultilevel, you need to remove the nperm argument in the fgsea function
## call.
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (5.94% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
## Warning in fgseaSimple(...): There were 14 pathways for which P-values were not
## calculated properly due to unbalanced gene-level statistic values
eaRes = eaRes %>% arrange(pval)

head(eaRes[,c(1:3,8)] , 40)
##                                                                                   pathway
##  1:                                            GOBP_VESICLE_MEDIATED_TRANSPORT_IN_SYNAPSE
##  2:                                                       GOBP_SYNAPTIC_VESICLE_RECYCLING
##  3:                                                          GOBP_PRESYNAPTIC_ENDOCYTOSIS
##  4:                                         GOBP_REGULATION_OF_SYNAPTIC_VESICLE_RECYCLING
##  5:                                       GOBP_REGULATION_OF_SYNAPTIC_VESICLE_ENDOCYTOSIS
##  6:                                                          GOBP_PHOSPHATE_ION_TRANSPORT
##  7:                                               GOBP_FATTY_ACID_TRANSMEMBRANE_TRANSPORT
##  8:                                              GOBP_L_GLUTAMATE_TRANSMEMBRANE_TRANSPORT
##  9:                                           GOBP_POSITIVE_REGULATION_OF_RRNA_PROCESSING
## 10:                                                      GOBP_ACIDIC_AMINO_ACID_TRANSPORT
## 11:                                                              GOBP_GLUTAMATE_SECRETION
## 12:                                       GOBP_L_ALPHA_AMINO_ACID_TRANSMEMBRANE_TRANSPORT
## 13:                                                           GOBP_L_AMINO_ACID_TRANSPORT
## 14:                                               GOBP_AMINO_ACID_TRANSMEMBRANE_TRANSPORT
## 15:                                                             GOBP_AMINO_ACID_TRANSPORT
## 16:                                             GOBP_ORGANIC_ACID_TRANSMEMBRANE_TRANSPORT
## 17:  GOBP_POSITIVE_REGULATION_OF_TRANSCRIPTION_ELONGATION_FROM_RNA_POLYMERASE_II_PROMOTER
## 18:                                                    GOBP_REGULATION_OF_RRNA_PROCESSING
## 19:                                                     GOBP_CELLULAR_RESPONSE_TO_ETHANOL
## 20:                                                      GOBP_DICARBOXYLIC_ACID_TRANSPORT
## 21:                                    GOBP_MAMMARY_GLAND_EPITHELIAL_CELL_DIFFERENTIATION
## 22:           GOBP_ANTIGEN_PROCESSING_AND_PRESENTATION_OF_PEPTIDE_ANTIGEN_VIA_MHC_CLASS_I
## 23: GOBP_ANTIGEN_PROCESSING_AND_PRESENTATION_OF_EXOGENOUS_PEPTIDE_ANTIGEN_VIA_MHC_CLASS_I
## 24:                                                           GOBP_CELL_REDOX_HOMEOSTASIS
## 25:                         GOBP_POSITIVE_REGULATION_OF_EXCITATORY_POSTSYNAPTIC_POTENTIAL
## 26:                                                       GOBP_NEUROTRANSMITTER_TRANSPORT
## 27:                                               GOBP_SYNAPTIC_TRANSMISSION_DOPAMINERGIC
## 28:                                       GOBP_C_TERMINAL_PROTEIN_AMINO_ACID_MODIFICATION
## 29:                                                            GOBP_UTP_METABOLIC_PROCESS
## 30:                           GOBP_G_PROTEIN_COUPLED_GLUTAMATE_RECEPTOR_SIGNALING_PATHWAY
## 31:                                                                GOBP_RESPIRATORY_BURST
## 32:       GOBP_NEGATIVE_REGULATION_OF_EPIDERMAL_GROWTH_FACTOR_ACTIVATED_RECEPTOR_ACTIVITY
## 33:                                            GOBP_NEUROTRANSMITTER_BIOSYNTHETIC_PROCESS
## 34:                                                 GOBP_CALCIUM_ION_REGULATED_EXOCYTOSIS
## 35:                                                     GOBP_ECTODERMAL_PLACODE_FORMATION
## 36:                                  GOBP_MODULATION_OF_EXCITATORY_POSTSYNAPTIC_POTENTIAL
## 37:                                                        GOBP_INORGANIC_ANION_TRANSPORT
## 38:                                                           GOBP_CAMP_METABOLIC_PROCESS
## 39:                                                      GOBP_PROTEIN_AUTOPHOSPHORYLATION
## 40:                                              GOBP_SYNAPTIC_TRANSMISSION_GLUTAMATERGIC
##                                                                                   pathway
##             pval      padj                               leadingEdge
##  1: 0.0004999900 0.8412843                            Slc17a7,Rab11a
##  2: 0.0005399892 0.8412843                              Slc17a7,Sncg
##  3: 0.0008199836 0.8412843                              Slc17a7,Sncg
##  4: 0.0008608781 0.8412843                             Slc17a7,Lrrk2
##  5: 0.0008681257 0.8412843                             Slc17a7,Lrrk2
##  6: 0.0020611955 0.8412843                           Slc17a7,Slc17a1
##  7: 0.0023400000 0.8412843                           Slc17a7,Slc17a8
##  8: 0.0029441807 0.8412843                           Slc17a7,Slc17a8
##  9: 0.0032649254 0.8412843 Trmt112,Wdr43,Riok2,Dimt1,Utp15,Riok1,...
## 10: 0.0041799164 0.8412843                            Slc17a7,Adora1
## 11: 0.0042602556 0.8412843                            Slc17a7,Adora1
## 12: 0.0046399072 0.8412843                              Slc17a7,Ace2
## 13: 0.0049199016 0.8412843                              Slc17a7,Ace2
## 14: 0.0053398932 0.8412843                              Slc17a7,Ace2
## 15: 0.0054798904 0.8412843                              Slc17a7,Ace2
## 16: 0.0057798844 0.8412843                             Slc17a7,Folr1
## 17: 0.0062421973 0.8412843   Tcerg1,Gtf2f1,Brd4,Eapp,Cdk13,Cdk12,...
## 18: 0.0064102564 0.8412843 Trmt112,Wdr43,Kat2b,Riok2,Dimt1,Utp15,...
## 19: 0.0064401150 0.8412843                               Cybb,Kcnmb1
## 20: 0.0077198456 0.8412843                             Slc17a7,Folr1
## 21: 0.0091099498 0.8412843                                Irf6,Erbb4
## 22: 0.0117797644 0.8412843                                 Cybb,Ncf2
## 23: 0.0121797564 0.8412843                                 Cybb,Ncf2
## 24: 0.0140800000 0.8412843                                 Cybb,Ncf2
## 25: 0.0152873379 0.8412843                              Chrna7,Rims1
## 26: 0.0152996940 0.8412843                 Slc17a7,Slc6a4,Grm4,Gpm6b
## 27: 0.0153159672 0.8412843                             Slc6a4,Chrnb2
## 28: 0.0153367324 0.8412843                               Agbl1,Irgm2
## 29: 0.0153917910 0.8412843        Entpd7,Ak3,Nme5,Nme4,Nme3,Nme1,...
## 30: 0.0154338843 0.8412843                                 Grm4,Grm8
## 31: 0.0154846454 0.8412843                                 Cybb,Ncf2
## 32: 0.0157155008 0.8412843                                Cblc,Zgpat
## 33: 0.0161157025 0.8412843                            Slc6a4,Slc44a4
## 34: 0.0167596648 0.8412843                               Scin,Notch1
## 35: 0.0187011576 0.8412843                                 Nrg3,Tbx3
## 36: 0.0188800000 0.8412843                              Chrna7,Rims1
## 37: 0.0194196116 0.8412843                              Slc17a7,Prnp
## 38: 0.0200084119 0.8412843                               Pde4c,Adcy4
## 39: 0.0201795964 0.8412843                         Tnk1,Epha1,Pdgfra
## 40: 0.0207595848 0.8412843                        Slc17a7,Grm4,Ntrk1
##             pval      padj                               leadingEdge

\(~\)
\(~\)

Differntial expression between conditions for each gene

Exploratory data visualization for some Nephron markers.

plotSmoothers(GAM.1to18, assays(GAM.1to18)$counts, gene = "Six2", alpha = 1, border = TRUE) + ggtitle("Six2")
grid.text("(Fig38)", 0.15, 0.85)

plotSmoothers(GAM.1to18, assays(GAM.1to18)$counts, gene = "Cited1", alpha = 1, border = TRUE) + ggtitle("Cited1")
grid.text("(Fig39)", 0.15, 0.85)

plotSmoothers(GAM.1to18, assays(GAM.1to18)$counts, gene = "Crym", alpha = 1, border = TRUE) + ggtitle("Crym")
grid.text("(Fig40)", 0.15, 0.85)

plotSmoothers(GAM.1to18, assays(GAM.1to18)$counts, gene = "C1qtnf12", alpha = 1, border = TRUE) + ggtitle("C1qtnf12/Fam132a")
grid.text("(Fig41)", 0.15, 0.85)

plotSmoothers(GAM.1to18, assays(GAM.1to18)$counts, gene = "Wnt4", alpha = 1, border = TRUE) + ggtitle("Wnt4")
grid.text("(Fig42)", 0.15, 0.85)

plotSmoothers(GAM.1to18, assays(GAM.1to18)$counts, gene = "S100g", alpha = 1, border = TRUE) + ggtitle("S100g")
grid.text("(Fig43)", 0.15, 0.85)

plotSmoothers(GAM.1to18, assays(GAM.1to18)$counts, gene = "Umod", alpha = 1, border = TRUE) + ggtitle("Umod")
grid.text("(Fig44)", 0.15, 0.85)

\(~\)
\(~\)

To test differential expression between conditions, we use the “conditionTest” function implemented in tradeSeq. We performed the global (considering lineage1 and 2 together) and pairwise (each lineage separately) statistical tests and set the log fold change threshold at 1.2.

condRes <- conditionTest(GAM.1to18, global = T, pairwise = T , l2fc = log2(1.2))

# log2(1) = 0
# log2(1.2) = 0.2630344
# log2(1.4) = 0.4854268

\(~\)
\(~\)

condRes$padj <- p.adjust(condRes$pvalue, "BH")
condRes$padj_lineage1 <- p.adjust(condRes$pvalue_lineage1, "BH")
condRes$padj_lineage2 <- p.adjust(condRes$pvalue_lineage2, "BH")

mean(condRes$padj <= 0.05, na.rm = TRUE)
## [1] 0.0009775702
# DEGs between wild-type and Wnt11 conditions
DEGs_general = subset(condRes , padj < 0.05)
DEGs_general
##           waldStat df       pvalue waldStat_lineage1 df_lineage1
## Eif2s3x   36.28846  8 1.554845e-05          16.28014           4
## Xist    1319.19793 16 0.000000e+00         831.88335           8
## Pbdc1     93.78951 16 4.986012e-13          47.60188           8
## Hddc3    155.91543 10 0.000000e+00          82.94289           5
## Crebzf    89.11504  4 0.000000e+00          48.78038           2
## Gm35082  106.07208  6 0.000000e+00          51.02320           3
## Hbb-bt   181.70547 16 0.000000e+00         143.34267           8
## Hbb-bs   129.20132 16 0.000000e+00          47.78361           8
## Hbb-y    117.07229 16 0.000000e+00          16.15757           8
## Hba-x     65.32296 16 6.481468e-08          33.36763           8
## Hba-a1    83.99710 16 3.142919e-11          34.58898           8
## Hba-a2    83.24145 16 4.312095e-11          35.13281           8
## Col1a1    60.26111 16 4.729297e-07          30.96917           8
## Meg3      48.33821 16 4.202527e-05          24.81855           8
## Kdm5d     65.38066  6 3.606893e-12          37.19822           3
## Eif2s3y  154.58927  6 0.000000e+00          78.62146           3
## Uty       39.44500  4 5.637265e-08          26.00062           3
## Ddx3y    197.02849  4 0.000000e+00         103.88385           2
##         pvalue_lineage1 waldStat_lineage2 df_lineage2 pvalue_lineage2
## Eif2s3x    2.665397e-03          20.00832           4    4.975146e-04
## Xist       0.000000e+00         487.31458           8    0.000000e+00
## Pbdc1      1.177078e-07          46.18763           8    2.189526e-07
## Hddc3      2.220446e-16          72.97254           5    2.464695e-14
## Crebzf     2.555489e-11          40.33466           2    1.743571e-09
## Gm35082    4.836731e-11          55.47784           4    2.579947e-11
## Hbb-bt     0.000000e+00          38.36279           8    6.450895e-06
## Hbb-bs     1.086656e-07          81.41771           8    2.531308e-14
## Hbb-y      4.017956e-02         100.91472           8    0.000000e+00
## Hba-x      5.286772e-05          31.95533           8    9.487318e-05
## Hba-a1     3.176323e-05          49.40812           8    5.309282e-08
## Hba-a2     2.528935e-05          48.10864           8    9.418119e-08
## Col1a1     1.422894e-04          29.29193           8    2.817783e-04
## Meg3       1.668499e-03          23.51966           8    2.757364e-03
## Kdm5d      4.177800e-08          28.18244           3    3.325465e-06
## Eif2s3y    1.110223e-16          75.96781           3    2.220446e-16
## Uty        9.534539e-06          13.79223           2    1.011707e-03
## Ddx3y      0.000000e+00          93.14465           2    0.000000e+00
##                 padj padj_lineage1 padj_lineage2
## Eif2s3x 1.684080e-02  1.000000e+00  5.376962e-01
## Xist    0.000000e+00  0.000000e+00  0.000000e+00
## Pbdc1   9.180743e-10  2.156878e-04  3.657106e-04
## Hddc3   0.000000e+00  8.137491e-13  7.751288e-11
## Crebzf  0.000000e+00  7.804464e-08  4.004330e-06
## Gm35082 0.000000e+00  1.266118e-07  6.771625e-08
## Hbb-bt  0.000000e+00  0.000000e+00  8.465878e-03
## Hbb-bs  0.000000e+00  2.156878e-04  7.751288e-11
## Hbb-y   0.000000e+00  1.000000e+00  0.000000e+00
## Hba-x   7.956218e-05  6.458320e-02  1.162070e-01
## Hba-a1  4.822548e-08  4.477150e-02  1.083860e-04
## Hba-a2  6.107585e-08  3.861683e-02  1.730391e-04
## Col1a1  5.442534e-04  1.533712e-01  3.235696e-01
## Meg3    4.298952e-02  1.000000e+00  1.000000e+00
## Kdm5d   6.037610e-09  9.569250e-05  4.699905e-03
## Eif2s3y 0.000000e+00  5.085932e-13  1.019906e-12
## Uty     7.414212e-05  1.588281e-02  9.783209e-01
## Ddx3y   0.000000e+00  0.000000e+00  0.000000e+00
# DEGs between wild-type and Wnt11 conditions in lineage1
DEGs_lineage1 = subset(condRes , padj_lineage1 < 0.05)
DEGs_lineage1
##           waldStat df       pvalue waldStat_lineage1 df_lineage1
## Xist    1319.19793 16 0.000000e+00         831.88335           8
## Pbdc1     93.78951 16 4.986012e-13          47.60188           8
## Hddc3    155.91543 10 0.000000e+00          82.94289           5
## Crebzf    89.11504  4 0.000000e+00          48.78038           2
## Gm35082  106.07208  6 0.000000e+00          51.02320           3
## Hbb-bt   181.70547 16 0.000000e+00         143.34267           8
## Hbb-bs   129.20132 16 0.000000e+00          47.78361           8
## Hba-a1    83.99710 16 3.142919e-11          34.58898           8
## Hba-a2    83.24145 16 4.312095e-11          35.13281           8
## Kdm5d     65.38066  6 3.606893e-12          37.19822           3
## Eif2s3y  154.58927  6 0.000000e+00          78.62146           3
## Uty       39.44500  4 5.637265e-08          26.00062           3
## Ddx3y    197.02849  4 0.000000e+00         103.88385           2
##         pvalue_lineage1 waldStat_lineage2 df_lineage2 pvalue_lineage2
## Xist       0.000000e+00         487.31458           8    0.000000e+00
## Pbdc1      1.177078e-07          46.18763           8    2.189526e-07
## Hddc3      2.220446e-16          72.97254           5    2.464695e-14
## Crebzf     2.555489e-11          40.33466           2    1.743571e-09
## Gm35082    4.836731e-11          55.47784           4    2.579947e-11
## Hbb-bt     0.000000e+00          38.36279           8    6.450895e-06
## Hbb-bs     1.086656e-07          81.41771           8    2.531308e-14
## Hba-a1     3.176323e-05          49.40812           8    5.309282e-08
## Hba-a2     2.528935e-05          48.10864           8    9.418119e-08
## Kdm5d      4.177800e-08          28.18244           3    3.325465e-06
## Eif2s3y    1.110223e-16          75.96781           3    2.220446e-16
## Uty        9.534539e-06          13.79223           2    1.011707e-03
## Ddx3y      0.000000e+00          93.14465           2    0.000000e+00
##                 padj padj_lineage1 padj_lineage2
## Xist    0.000000e+00  0.000000e+00  0.000000e+00
## Pbdc1   9.180743e-10  2.156878e-04  3.657106e-04
## Hddc3   0.000000e+00  8.137491e-13  7.751288e-11
## Crebzf  0.000000e+00  7.804464e-08  4.004330e-06
## Gm35082 0.000000e+00  1.266118e-07  6.771625e-08
## Hbb-bt  0.000000e+00  0.000000e+00  8.465878e-03
## Hbb-bs  0.000000e+00  2.156878e-04  7.751288e-11
## Hba-a1  4.822548e-08  4.477150e-02  1.083860e-04
## Hba-a2  6.107585e-08  3.861683e-02  1.730391e-04
## Kdm5d   6.037610e-09  9.569250e-05  4.699905e-03
## Eif2s3y 0.000000e+00  5.085932e-13  1.019906e-12
## Uty     7.414212e-05  1.588281e-02  9.783209e-01
## Ddx3y   0.000000e+00  0.000000e+00  0.000000e+00
# DEGs between wild-type and Wnt11 conditions in lineage2
DEGs_lineage2 = subset(condRes , padj_lineage2 < 0.05)
DEGs_lineage2
##           waldStat df       pvalue waldStat_lineage1 df_lineage1
## Xist    1319.19793 16 0.000000e+00      831.88335345           8
## Pbdc1     93.78951 16 4.986012e-13       47.60188370           8
## Hddc3    155.91543 10 0.000000e+00       82.94288767           5
## Crebzf    89.11504  4 0.000000e+00       48.78038163           2
## Gm35082  106.07208  6 0.000000e+00       51.02319583           3
## Hbb-bt   181.70547 16 0.000000e+00      143.34267459           8
## Hbb-bs   129.20132 16 0.000000e+00       47.78361213           8
## Hbb-y    117.07229 16 0.000000e+00       16.15756801           8
## Nktr      28.70845  6 6.907214e-05        0.04676395           3
## Hba-a1    83.99710 16 3.142919e-11       34.58897703           8
## Hba-a2    83.24145 16 4.312095e-11       35.13280800           8
## Kdm5d     65.38066  6 3.606893e-12       37.19821884           3
## Eif2s3y  154.58927  6 0.000000e+00       78.62145884           3
## Ddx3y    197.02849  4 0.000000e+00      103.88384644           2
##         pvalue_lineage1 waldStat_lineage2 df_lineage2 pvalue_lineage2
## Xist       0.000000e+00         487.31458           8    0.000000e+00
## Pbdc1      1.177078e-07          46.18763           8    2.189526e-07
## Hddc3      2.220446e-16          72.97254           5    2.464695e-14
## Crebzf     2.555489e-11          40.33466           2    1.743571e-09
## Gm35082    4.836731e-11          55.47784           4    2.579947e-11
## Hbb-bt     0.000000e+00          38.36279           8    6.450895e-06
## Hbb-bs     1.086656e-07          81.41771           8    2.531308e-14
## Hbb-y      4.017956e-02         100.91472           8    0.000000e+00
## Nktr       9.973478e-01          28.66169           3    2.637626e-06
## Hba-a1     3.176323e-05          49.40812           8    5.309282e-08
## Hba-a2     2.528935e-05          48.10864           8    9.418119e-08
## Kdm5d      4.177800e-08          28.18244           3    3.325465e-06
## Eif2s3y    1.110223e-16          75.96781           3    2.220446e-16
## Ddx3y      0.000000e+00          93.14465           2    0.000000e+00
##                 padj padj_lineage1 padj_lineage2
## Xist    0.000000e+00  0.000000e+00  0.000000e+00
## Pbdc1   9.180743e-10  2.156878e-04  3.657106e-04
## Hddc3   0.000000e+00  8.137491e-13  7.751288e-11
## Crebzf  0.000000e+00  7.804464e-08  4.004330e-06
## Gm35082 0.000000e+00  1.266118e-07  6.771625e-08
## Hbb-bt  0.000000e+00  0.000000e+00  8.465878e-03
## Hbb-bs  0.000000e+00  2.156878e-04  7.751288e-11
## Hbb-y   0.000000e+00  1.000000e+00  0.000000e+00
## Nktr    6.693817e-02  1.000000e+00  4.038426e-03
## Hba-a1  4.822548e-08  4.477150e-02  1.083860e-04
## Hba-a2  6.107585e-08  3.861683e-02  1.730391e-04
## Kdm5d   6.037610e-09  9.569250e-05  4.699905e-03
## Eif2s3y 0.000000e+00  5.085932e-13  1.019906e-12
## Ddx3y   0.000000e+00  0.000000e+00  0.000000e+00
# Common DEGs between the two lineages
intersect(DEGs_lineage1 , DEGs_lineage2)
##           waldStat df       pvalue waldStat_lineage1 df_lineage1
## Xist    1319.19793 16 0.000000e+00         831.88335           8
## Pbdc1     93.78951 16 4.986012e-13          47.60188           8
## Hddc3    155.91543 10 0.000000e+00          82.94289           5
## Crebzf    89.11504  4 0.000000e+00          48.78038           2
## Gm35082  106.07208  6 0.000000e+00          51.02320           3
## Hbb-bt   181.70547 16 0.000000e+00         143.34267           8
## Hbb-bs   129.20132 16 0.000000e+00          47.78361           8
## Hba-a1    83.99710 16 3.142919e-11          34.58898           8
## Hba-a2    83.24145 16 4.312095e-11          35.13281           8
## Kdm5d     65.38066  6 3.606893e-12          37.19822           3
## Eif2s3y  154.58927  6 0.000000e+00          78.62146           3
## Ddx3y    197.02849  4 0.000000e+00         103.88385           2
##         pvalue_lineage1 waldStat_lineage2 df_lineage2 pvalue_lineage2
## Xist       0.000000e+00         487.31458           8    0.000000e+00
## Pbdc1      1.177078e-07          46.18763           8    2.189526e-07
## Hddc3      2.220446e-16          72.97254           5    2.464695e-14
## Crebzf     2.555489e-11          40.33466           2    1.743571e-09
## Gm35082    4.836731e-11          55.47784           4    2.579947e-11
## Hbb-bt     0.000000e+00          38.36279           8    6.450895e-06
## Hbb-bs     1.086656e-07          81.41771           8    2.531308e-14
## Hba-a1     3.176323e-05          49.40812           8    5.309282e-08
## Hba-a2     2.528935e-05          48.10864           8    9.418119e-08
## Kdm5d      4.177800e-08          28.18244           3    3.325465e-06
## Eif2s3y    1.110223e-16          75.96781           3    2.220446e-16
## Ddx3y      0.000000e+00          93.14465           2    0.000000e+00
##                 padj padj_lineage1 padj_lineage2
## Xist    0.000000e+00  0.000000e+00  0.000000e+00
## Pbdc1   9.180743e-10  2.156878e-04  3.657106e-04
## Hddc3   0.000000e+00  8.137491e-13  7.751288e-11
## Crebzf  0.000000e+00  7.804464e-08  4.004330e-06
## Gm35082 0.000000e+00  1.266118e-07  6.771625e-08
## Hbb-bt  0.000000e+00  0.000000e+00  8.465878e-03
## Hbb-bs  0.000000e+00  2.156878e-04  7.751288e-11
## Hba-a1  4.822548e-08  4.477150e-02  1.083860e-04
## Hba-a2  6.107585e-08  3.861683e-02  1.730391e-04
## Kdm5d   6.037610e-09  9.569250e-05  4.699905e-03
## Eif2s3y 0.000000e+00  5.085932e-13  1.019906e-12
## Ddx3y   0.000000e+00  0.000000e+00  0.000000e+00

\(~\)
\(~\)
\(~\)
\(~\)

Bootstrap Hypothesis Test

First, we performed a t-test and Wilcoxon/Mann-whitney test on the distribution of pseudotime1 between wild-type and wnt11 groups. Both of which resulted in a p-value less than 0.05 percent.

# Comparing the means of two distribution with t test
t.test(slingPseudotime(slingshot.wild.1to18)[, 1], slingPseudotime(slingshot.wnt11.1to18)[, 1])
## 
##  Welch Two Sample t-test
## 
## data:  slingPseudotime(slingshot.wild.1to18)[, 1] and slingPseudotime(slingshot.wnt11.1to18)[, 1]
## t = -5.1732, df = 1721.2, p-value = 2.57e-07
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.9289447 -0.4181948
## sample estimates:
## mean of x mean of y 
##  3.274941  3.948511
# Wilcoxon/Mann-whitney compares the medians of two distributions
wilcox.test(slingPseudotime(slingshot.wild.1to18)[, 1], slingPseudotime(slingshot.wnt11.1to18)[, 1])
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  slingPseudotime(slingshot.wild.1to18)[, 1] and slingPseudotime(slingshot.wnt11.1to18)[, 1]
## W = 322376, p-value = 1.517e-06
## alternative hypothesis: true location shift is not equal to 0

\(~\)
\(~\)

Performing two Bootstrap hypothesis test based on mean and median of the distributions.

# In Bootstrap, we consider the null hypothesis is true, that is , the probability that
# a specific pseudotime comes from wild type is the same as it comes from Wnt11. Therefore,
# resampling is applied on aggregated data.

wild_pt1 = as.numeric(na.omit(slingPseudotime(slingshot.wild.1to18)[, 1]))
wnt11_pt1 = as.numeric(na.omit(slingPseudotime(slingshot.wnt11.1to18)[, 1]))

wild_pt1_mean = mean(wild_pt1)

wnt11_pt1_mean = mean(wnt11_pt1)

wild_pt1_median = median(wild_pt1)

wnt11_pt1_median = median(wnt11_pt1)

wild_pt1_mean
## [1] 3.274941
wnt11_pt1_mean
## [1] 3.948511
wild_pt1_median
## [1] 2.841504
wnt11_pt1_median
## [1] 3.375997
test_stat_mean = abs(wild_pt1_mean - wnt11_pt1_mean)
test_stat_median = abs(wild_pt1_median - wnt11_pt1_median)

L = length(c(wild_pt1, wnt11_pt1))
L
## [1] 1732
# Number of resampling 
B = 10000

numbers = c(wild_pt1, wnt11_pt1)

# resampling L*B times
Bootstrap_matrix = matrix(sample(x = numbers , size = L*B , replace = T) , nrow = L, ncol = B)
dim(Bootstrap_matrix)
## [1]  1732 10000
Boot_stat_mean = rep(0,B)
Boot_stat_median = rep(0,B)

length(wild_pt1)
## [1] 792
for(i in 1:B){
        
        Boot_stat_mean[i] = abs(mean(Bootstrap_matrix[1:792, i]) - mean(Bootstrap_matrix[793:1732, i]))
        Boot_stat_median[i] = abs(median(Bootstrap_matrix[1:792, i]) - median(Bootstrap_matrix[793:1732, i]))
        
}

# Density plot of Bootstrap statistics
Boot_stat_mean[1:20]
##  [1] 0.08763163 0.08518836 0.15115668 0.03502281 0.05215037 0.14190328
##  [7] 0.17647074 0.02751870 0.06701882 0.08203811 0.13314141 0.05520239
## [13] 0.01391266 0.20988616 0.22387382 0.02746928 0.13939191 0.05610408
## [19] 0.18180438 0.03027709
dens = density(Boot_stat_mean)
plot(dens, frame =  F , col = "steelblue", main = "Bootstrap statistics means") 
polygon(dens, col = "steelblue")
grid.text("(Fig45)", 0.15, 0.75)

Boot_stat_median[1:20]
##  [1] 0.03267741 0.32980174 0.25902104 0.12543396 0.09649757 0.03764550
##  [7] 0.50266934 0.06891447 0.11483682 0.25922714 0.24215826 0.12469570
## [13] 0.04605807 0.20465785 0.23964808 0.04441800 0.12099561 0.18582881
## [19] 0.24390804 0.09723747
dens = density(Boot_stat_median)
plot(dens, frame =  F , col = "indianred1", main = "Bootstrap statistics medians") 
polygon(dens, col = "indianred1")
grid.text("(Fig46)", 0.15, 0.75)

# Computing the p-value for mean statistics and median statistics. 
# Both Hypothesis tests gave rise to a p-value less than 0.05 percent. 

sum(Boot_stat_mean > test_stat_mean)
## [1] 0
sum(Boot_stat_median > test_stat_median)/10000
## [1] 0.0329
mean(Boot_stat_median > test_stat_median)
## [1] 0.0329

\(~\)
\(~\)

Differential Expression Between Clusters

Differential expression between Wnt11KO vs Wild-type in all clusters. 1_wnt11-hom: wnt11 cells in cluster1 and so on.

so.small$Cluster_Genotype = paste(Idents(so.small), so.small$genotype, sep = "_")
Idents(so.small) = "Cluster_Genotype"

DEGs = FindMarkers(so.small, ident.1 = "1_wnt11-hom", ident.2 = "1_wild-type", verbose = F)
DEGs %>% arrange(desc(abs(avg_log2FC))) %>% filter(abs(avg_log2FC) > 0.5 & p_val_adj < 0.05)
##                  p_val avg_log2FC pct.1 pct.2     p_val_adj
## Xist     1.178851e-162 -3.3946361 0.056 1.000 2.186650e-158
## Gm35082   2.958657e-58  0.7788606 0.458 0.012  5.488013e-54
## Ddx3y     6.301919e-61  0.7426733 0.500 0.027  1.168943e-56
## Eif2s3y   9.122415e-53  0.6252948 0.432 0.017  1.692117e-48
## Ftl1-ps1  1.113321e-10 -0.5785791 0.596 0.762  2.065100e-06
DEGs = FindMarkers(so.small, ident.1 = "21_wnt11-hom", ident.2 = "21_wild-type", verbose = F)
DEGs %>% arrange(desc(abs(avg_log2FC))) %>% filter(abs(avg_log2FC) > 0.5 & p_val_adj < 0.05)
##             p_val avg_log2FC pct.1 pct.2    p_val_adj
## Xist 1.019964e-09  -3.990642     0     1 1.891931e-05
DEGs = FindMarkers(so.small, ident.1 = "5_wnt11-hom", ident.2 = "5_wild-type", verbose = F)
DEGs %>% arrange(desc(abs(avg_log2FC))) %>% filter(abs(avg_log2FC) > 0.5 & p_val_adj < 0.05)
##                 p_val avg_log2FC pct.1 pct.2    p_val_adj
## Xist    2.186199e-100 -3.8232985 0.007 0.995 4.055181e-96
## Ddx3y    5.854359e-37  0.8172825 0.586 0.010 1.085925e-32
## Crabp1   2.507432e-08  0.8153917 0.500 0.285 4.651036e-04
## Cdkn1c   1.169589e-06  0.6940403 0.967 0.938 2.169471e-02
## Eif2s3y  3.342528e-29  0.6736892 0.516 0.026 6.200055e-25
## Gm35082  1.012469e-28  0.6589631 0.467 0.000 1.878028e-24
## Meg3     1.493468e-06  0.5397005 0.556 0.326 2.770234e-02
## Hddc3    2.344082e-15  0.5089163 0.539 0.192 4.348038e-11
#subset(DEGs , abs(avg_log2FC) > 0.5 & p_val_adj < 0.05)

DEGs = FindMarkers(so.small, ident.1 = "13_wnt11-hom", ident.2 = "13_wild-type", verbose = F)
DEGs %>% arrange(desc(abs(avg_log2FC))) %>% filter(abs(avg_log2FC) > 0.5 & p_val_adj < 0.05)
##                 p_val avg_log2FC pct.1 pct.2    p_val_adj
## Xist     1.586450e-38 -3.1064394 0.033 0.986 2.942706e-34
## Ftl1-ps1 1.476746e-13 -1.6014272 0.388 0.838 2.739216e-09
## Ddx3y    4.219953e-12  0.7150880 0.512 0.027 7.827591e-08
## Gm35082  7.515481e-11  0.6509924 0.430 0.000 1.394047e-06
## Eif2s3y  5.504918e-10  0.5853334 0.421 0.014 1.021107e-05
DEGs = FindMarkers(so.small, ident.1 = "18_wnt11-hom", ident.2 = "18_wild-type", verbose = F)
DEGs %>% arrange(desc(abs(avg_log2FC))) %>% filter(abs(avg_log2FC) > 0.5 & p_val_adj < 0.05)
##             p_val avg_log2FC pct.1 pct.2    p_val_adj
## Xist 9.585284e-18  -3.058509 0.043     1 1.777974e-13
DEGs = FindMarkers(so.small, ident.1 = "15_wnt11-hom", ident.2 = "15_wild-type", verbose = F)
DEGs %>% arrange(desc(abs(avg_log2FC))) %>% filter(abs(avg_log2FC) > 0.5 & p_val_adj < 0.05)
##             p_val avg_log2FC pct.1 pct.2    p_val_adj
## Xist 9.231529e-23  -2.971087  0.05     1 1.712356e-18